Probabilistic bounded reachability for hybrid systems with continuous nondeterministic and probabilistic parameters
Fedor Shmarov
Paolo Zuliani
%Verification of cyber-physical systems is a difficult, yet extremely important, problem. %Hybrid systems offer a modelling framework in which to perform formal verification of %cyber-physical systems. We develop an algorithm for computing bounded reachability probability for hybrid systems, i.e., the probability that the system reaches an unsafe region within a finite number of discrete transitions. In particular, we focus on hybrid systems with continuous dynamics given by solutions of nonlinear ordinary differential equations (with possibly nondeterministic initial conditions and parameters), and probabilistic behaviour given by initial parameters distributed as continuous (with possibly infinite support) and discrete random variables. Our approach is to define an appropriate relaxation of the (undecidable) reachability problem, so that it can be solved by $\delta$-complete decision procedures. In particular, for systems with continuous random parameters only, we develop a validated integration procedure which computes an arbitrarily small interval that is guaranteed to contain the reachability probability. In the more general case of systems with both nondeterministic and probabilistic parameters, our procedure computes a guaranteed enclosure for the range of reachability probabilities. We have applied our approach to a number of nonlinear hybrid models and validated the results by comparison with Monte Carlo simulation. %Moreover, we discuss certain types of random variables with unbounded support and show %that the bounded $\delta$-reachability problem can still be solved by using an appropriate %$\delta$-complete decision procedure.
algorithm.tex

\section{Computing $\delta$-Reachability Probability}

\subsection{Computing indicator functions}
From Algorithm \ref{alg:integration} we obtain a partition of the domain of the random
parameters which will guarantee the computation of integral (\ref{eq:probabilistic-reachability_rv})
with the desidered accuracy. In general, given $z_0\in Z$, the reachability probability is computed 
by integrating the probability measure of the random parameters over the restriction $B|_{z_0}$. 
We need to compute the following integral
\[
\int_{B|_{z_0}} \,dP(r)\quad  \left(= \int_{D\times R} 
        I_{U}(r, z_0)dP(r)\right)
\]
where $B$ is the set (\ref{def:Borel-set}), $z_0\in Z$, and $I_{U}$ is the indicator function
\begin{equation*}
        I_{U}(r, z_0) = 
        \begin{cases}
                1 \quad \text{if the system with parameter $(r, z_0)$ reaches
                        $U$ in $k$ steps} \\
                0  \quad\text{otherwise.}
        \end{cases}
\end{equation*}
We now show how to compute $I_U$ or, equivalently, set $B$.
Let $[\rho]\subseteq \Lambda$ be a box and $\phi$ be a formula of the form:
\begin{equation} \label{eq:phi}
\phi([\rho]) = \exists \lambda \in [\rho]\ \psi(\lambda)
\end{equation}
If the formula is true then $[\rho]$ contains a value for the initial parameters for which
the system reaches the unsafe region $U$.
Taking the complement of the unsafe region $U^{C} = S / U$ ($S$ is the state space of the 
system) and defining a predicate $\text{unsafe}^{C}_{q}(\textbf{x}^{t}) \equiv 
((q, \textbf{x}^{t}_{q}) \in U^{C})$ we want to ensure that the system never reaches the 
unsafe region {\em within} the $k$-th step with an initial parameter from $[\rho]$. 
In order to conclude that it is sufficient to evaluate the formula:
\begin{equation} \label{eq:phi_C}
\begin{split}
\phi^{C}([\rho]) = 
	\exists \lambda \in [\rho], \exists \textbf{x}^{0}_{0,q_0}, \exists \textbf{x}^{t}_{0,q_0}, \exists t_{0, q_0} , ... , \exists \textbf{x}^{0}_{0,q_m}, \\
\exists \textbf{x}^{t}_{0,q_m}, \exists t_{0, q_m}, ... , \exists \textbf{x}^{0}_{k,q_m}, \exists \textbf{x}^t_{k,q_m}, \exists t_{k, q_m}, \forall t'_{k, q_m} \in [0, t_{k, q_m}]:\\
(\bigvee_{q \in Q} (\text{init}_{q}(\textbf{x}^{0}_{0,q}) \wedge \text{flow}_{q}(\lambda, \textbf{x}^{0}_{0,q}, \textbf{x}^{t}_{0,q}))) \wedge \\
(\bigwedge_{i=0}^{k-1} (\bigvee_{q,q' \in Q}(\text{jump}_{q \rightarrow q'}(\lambda, \textbf{x}^{t}_{i,q}, \textbf{x}^{0}_{i+1,q'}) 
\wedge \\
(\text{flow}_{q'}(\lambda, \textbf{x}^{0}_{i+1,q'}, \textbf{x}^{t}_{i+1,q'}))) \wedge \\
(\bigvee_{q \in Q} (\text{unsafe}^{C}_{q}(\textbf{x}^{t}_{k,q}) \wedge (\text{jump}_{q \rightarrow q'}(\lambda, \textbf{x}^{t}_{k,q}, \textbf{x}^{0}_{k+1,q'}) \vee 
(t_{k, q_m} \ge T))))))
\end{split}
\end{equation}
Note that $\phi^C$ is not the logical negation of $\phi$ --- it is in fact an $\exists\forall$-quantified
formula.
The last term of $\phi^C$ ensures that the system either does not reach the unsafe region on 
the $k$-th step before it can make a transition to the successor mode or it reaches the time bound before reaching the unsafe region. This should not be confused with reaching the time bound in any of the preceding modes as it means that the system fails to reach the $k$-th step and should be, therefore, unsatisfiable.
If the formula evaluates to true then the system does not reach the unsafe region on the $k$-th step.
Then, set $B$ can be defined as a finite collection
$\{[\rho]_{i} : \phi([\rho]_{i}) \wedge (\neg \phi^{C}([\rho]_{i})\}$.
\hide{
By integrating a probability density function of a random variable we obtain intervals 
such that size of the enclosure containing the exact value of a partial sum for each 
of them is smaller than a local error.
Then a decision about the relation between each of the obtained intervals and the (unknown) 
Borel set $B$ should be taken. In order to achieve this we use a $\delta$-complete decision 
procedure, \ie, dReal \cite{DBLP:conf/cade/GaoKC13}, to verify formulas (\ref{eq:phi}) and 
(\ref{eq:phi_C}). 
\iffalse
The outcomes of the formula evaluations are shown in the diagram (Fig. \ref{fig:diagram}):

\begin{figure} [ht!]
\centering 
\includegraphics[width=70mm]{diagram.pdf}
\caption{The analysis of the interval $[r]$ from the domain of random variable}
\end{figure} \label{fig:diagram}
\fi
}
To build such a collection, we iteratively evaluate $\phi$ and $\phi^C$ with a $\delta$-complete
procedure (\eg, dReal \cite{DBLP:conf/cade/GaoKC13}). Given a box $[\rho]$, there are four 
possible outcomes:
\begin{itemize}
	\item{$\phi([\rho])$ is \textbf{unsat}. Hence, there are {\em for sure} no values in $[\rho]$ 
	such that the system reaches the unsafe region, so $[\rho]$ is not in $B$.}
	\item{$\phi([\rho])$ is \textbf{$\delta$-sat}. Then, there is a value in $[\rho]$ such that 
	the system reaches $U$ or $U^{\delta}$ ($\delta$-weakening of set $U$).}
	\item{$\phi^C([\rho])$ is \textbf{unsat}. Therefore, there is {\em for sure} no value in 
	$[\rho]$ such that for all time points on the $k$-th step the system stays in $U^C$. 
	In other words, for all the values in $[\rho]$ the system 
	reaches $U$, so $[\rho]$ is fully contained in $B$.}
	\item{$\phi^C([\rho])$ is \textbf{$\delta$-sat}. Then there is a value in $[\rho]$ such that the system stays within $U^{C}$ or $U^{C^{\delta}}$. In combination with outcome \textbf{$\delta$-sat} for $\phi([\rho])$ it signals that $[\rho]$ is a {\em mixed} interval (it contains values from both 
$B$ and $B^C$).}
\end{itemize}
Therefore, \textbf{unsat} answers enable us to decide whether $[\rho]$ is a subset of or disjoint from
set $B$. If \textbf{$\delta$-sat} is returned for both formulae, then we are either dealing with a 
false alarm (an unsatisfiable formula is verified as \textbf{$\delta$-sat} because of the
overapproximation) or a mixed interval.

\vspace{-2ex}
\subsection{Main algorithm}
\vspace{-1ex}
The overapproximation (controlled by $\delta$) introduced by $\delta$-complete procedures can 
cause false alarms. We thus begin by addressing the choice of $\delta$.
Obviously, it is impossible to decide {\em correctly} (\ie, obtaining \textbf{unsat} for one of 
$\phi$ and $\phi^C$) on each interval if a fixed $\delta$ (even a very small one) is used for
evaluate all formulae. 

%For example, consider the formula $\phi : \exists x \in [0, \frac{\pi}{2}]: (\cos{x}-0.5=0)$ and let $\delta=1$. Then a $\delta$-weakening of presented formula is $\phi^{\delta} : \exists x \in [0, \frac{\pi}{2}]: |\cos{x}-0.5| \le 1$.
%$\phi$ is true on the interval $[\frac{\pi}{4}, \frac{\pi}{2}]$ and false on $[0, \frac{\pi}{4}]$. However, $\phi^{\delta}$ is true for all the values on the interval $[0, \frac{\pi}{2}]$ and it is not possible to decide correctly on a smaller interval with the same or greater value of $\delta$. However, reducing $\delta$ (for instance, twice for each next iteration) there will be obtained intervals where the formula is \textbf{unsat}.

\begin{lemma} \label{lemma:converge}
Let $\phi$ be an arbitrary bounded $\Sigma_1$ formula and $\phi^{\delta}$ its weakening. Then the following holds:
\begin{equation*}
\forall \delta, \delta' \in \mathbb{Q}^{+}, 0 \leq \delta' < \delta: \neg \phi^{\delta} \rightarrow \neg \phi^{\delta'} \rightarrow  \neg \phi
\end{equation*}
\end{lemma}
(See Appendix \ref{apndx:delta-sat} for an overview of $\delta$-weakening.)
Lemma \ref{lemma:converge} means that unsatisfiability of a weakened formula 
implies unsatisfiability of its strengthening and of the initial formula. 
We next show that when an interval is uncertain, by applying Lemma \ref{lemma:converge} 
we can obtain $\delta$ and a subinterval for which 
a $\delta$-complete decision procedure can give a correct answer.
\begin{proposition} \label{prop:converge}
Let $\phi$ and $\phi^C$ as per (\ref{eq:phi}) and (\ref{eq:phi_C}), and $[u,v]$ an interval. Then:
\begin{equation*}
\begin{split}
&\exists \delta \in \mathbb{Q}^{+}:(\phi([u, v])-\delta\textbf{-sat}) \wedge (\phi^{C}([u, v])-\delta\textbf{-sat}) \Rightarrow \\
&\exists [u', v'] \subseteq [u, v]: (\phi([u', v'])-\textbf{unsat}) \oplus (\phi^{C}([u', v'])-\textbf{unsat})
\end{split}
\end{equation*}
where $\oplus$ denotes exclusive or.
\end{proposition}

%\begin{proposition}\label{prop:bounded}
%Given $\epsilon \in (0, 1] \cap \mathbb{Q}$, $k\in\mathbb{N}$ and a hybrid system with one 
%\textbf{bounded} continuous random initial parameter, there exists an algorithm for computing an interval 
%of size not larger than $\epsilon$ that contains the value of (\ref{eq:probabilistic-reachability}), 
%\ie, the probability of reaching the unsafe region in $k$ steps.
%\end{proposition}

We now present the full algorithm for computing bounded reachability probability. We begin by addressing
random initial parameters with (un)bounded support. Given $\epsilon \in (0,1]\cap\mathbb{Q}$, it is always 
to possible to find a {\em bounded} region of the random variable support with area larger
than $1-\epsilon$. In fact, such a problem can be stated as a $\delta$-satisfiability question 
and thus solved by a $\delta$-complete procedure. Therefore, the verified integration procedure presented
in Section \ref{sec:vip} can be applied to a random variable with unbounded domain.
If we introduce multiple independent random parameters we can still use the same verified 
integration procedure provided that each random variable is integrated with a higher accuracy, as the
next proposition shows.
\begin{proposition}\label{prop:multrv}
%To compute the integral of $l$ independent continuous random initial parameters within an interval of 
%length $\epsilon_{prod} \in (0,1]\cap\mathbb{Q}$, it is sufficient that each random variable is integrated 
%with precision $\epsilon$ satisfying:
Given a hybrid system with $l$ independent continuous random parameters, to compute with precision 
$\epsilon_{prod} \in (0,1]\cap\mathbb{Q}$ the reachability probability it is sufficient that each random 
variable is integrated with precision $\epsilon$ satisfying:
\begin{equation} \label{eq:multivarmaininequalityspecialcase}
\epsilon_{prod} \ge \sum_{i=1}^l {l\choose{i}}\epsilon^i 
\end{equation}
where ${l}\choose{i}$ is the binomial coefficient.
\end{proposition}

Suppose now a hybrid system has (continuous) nondeterministic parameters. Then the 
probability that the system reaches the unsafe region becomes a function of the nondeterministic 
parameters. In particular, the indicator function 
$I_{U}(r, z)$ can be equal to 0 and 1 for the same values of the continuous
{\em random} parameters, \ie, there may exist $r_0$ and $z_0 \neq z_1$ such that 
$I_{U}(r_0, z_0) = 0$ and $I_{U}(r_0, z_1) = 1$. 
Therefore, it is in general impossible to provide any guarantees on the length of probability interval,
and we need to compute an enclosure for all probabilities.
We will use the following symbolic notation for hybrid systems:
\begin{itemize}
	\item{\textbf{HA} ({\em Hybrid Automaton})} - a hybrid system without initial random parameters 
		(only deterministic and nondeterministic).
	\item{\textbf{PHA} ({\em Probabilistic Hybrid Automaton})} - a hybrid system with random and 
		deterministic continuous initial parameters (no nondeterminism).
	\item{\textbf{NPHA} ({\em Nondeterministic Probabilistic Hybrid Automaton})} - a hybrid system 
		with random, deterministic and nondeterministic continuous initial parameters.
\end{itemize}
We first state the algorithm for \textbf{NPHA}s with no discrete probability.
\begin{proposition}\label{prop:cn_nd}
Given $\epsilon \in (0, 1] \cap \mathbb{Q}$, $k\in\mathbb{N}$, and an {\em\textbf{NPHA}} {\em without discrete
random parameters}, there exists an algorithm for computing an interval containing the set of $k$-step 
reachability probabilities. If the system has no nondeterministic parameters, the algorithm returns an 
interval of size not larger than $\epsilon$ containing the $k$-step reachability probability 
(\ref{eq:probabilistic-reachability_rv}).
\end{proposition}
The pseudo-code of the algorithm is presented in Algorithm \ref{alg:prob-reach}. Informally,
the algorithm starts by getting an interval partition from the validated integration
procedure (Algorithm \ref{alg:integration}) for each random variable; also, a candidate
probability interval is initialised to $[0,1]$. 
Then, it evaluates the formulae $\phi$ and $\phi^C$ on the current partition, which will be 
refined whenever both $\phi$ and $\phi^C$ are $\delta$-\textbf{sat}. Instead, an \textbf{unsat} 
answer is used to refine the probability interval. The termination condition depends on the model type.
If there are no nondeterministic parameters, then the algorithm will terminate when the width 
of the probability interval satisfies the desired size $\epsilon$. Otherwise, the algorithm 
terminates when the maximum length of the boxes in the partition is smaller than $\epsilon$. 
(Given a box we can split it into $2^{n}$ boxes of (pairwise) equal size in such a way that each 
interval in the box is reduced. However, any division strategy can be 
applied as long as the size of each interval forming the box is reduced.)



\begin{algorithm}
\label{alg:prob-reach}
input: continuous random parameters $\overline{r} = \{r_{1}, ..., r_l\}$ with their probability 
densities $\overline{f}(r)$, $\epsilon_{prod}\in (0,1]\cap\mathbb{Q}$, hybrid system description $\phi$\;
output: interval $[I]$ enclosing the reachability probabilities\\
\tcp{obtain $\epsilon$ from (\ref{eq:multivarmaininequalityspecialcase}) using $\epsilon_{prod}$}
$\epsilon_{inf} = t \epsilon$\;
$\epsilon_{prob} = (1-t) \epsilon$\;

\tcp{obtain bounds from (\ref{eq:epsilon-inf}) for each continuous random parameter}
$[\overline{a}, \overline{b}] = \bigcup_{i=1}^{l}($\textbf{bounds}$(f(r_{i}), \epsilon_{inf})$)\;

$[\overline{r}] = \bigcup_{i=1}^{l}($\textbf{Algorithm1}$(f(r_{i}), [a_{i}, b_{i}], \epsilon_{prob})$)\;

\tcp{get Cartesian product of intervals from obtained partitions}
$B$.push$([r_{1}]\times \cdots \times[r_{l}])$\;
\tcp{set initial probability intervals for upper/lower approximation}
$[P_{lower}] = [0.0, 0.0]$\;
$[P_{upper}] = [1.0, 1.0]$\;
\tcp{consider unbounded segments}
$[P_{upper}] = [P_{upper}] + 1 - \prod_{i=1}^{l}\int_{a_{i}}^{b_{i}}\, f_{i}(x)\, dx $\;
\While{true}
{
	\tcp{stack containing extra divisions of the boxes}
	$D = \emptyset$\;
	\While{$size(B) > 0$}
	{
	$\textbf{box} = B$.pop()\;
	\tcp{$\delta$-complete procedure evaluates formula}
	\If{$\phi(\textbf{box})$ = $\delta$-sat}
	{
		\tcp{$\delta$-complete procedure evaluates formula}
		\If{$\phi^{C}(\textbf{box})$ = $\delta$-sat}
		{
			\tcp{split initial box into $2^{l}$ boxes of equal size}
			$D$.push(branch$(\textbf{box})$)\;
		}\Else
		{
			\tcp{increase lower approximation}
			$[P_{lower}] = [P_{lower}] + [S](\textbf{box})$\;
		}
	} \Else
	{
		\tcp{decrease upper approximation}
		$[P_{upper}] = [P_{upper}] - [S](\textbf{box})$\;
	}
	}
	$B = D$\;
	\If{MODEL\_TYPE($\phi$) = \textbf{PHA}}
	{
		\tcp{termination condition when nondeterminism is absent}
		\If{$\overline{[P_{upper}]} - \underline{[P_{lower}]} \le \epsilon_{prob}$}
		{
			\Return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]$\;
		}
	}
	\If{MODEL\_TYPE($\phi$) = \textbf{NPHA}}
	{
		\If{$\max_{\textbf{box} \in B}(|\textbf{box}|) \le \epsilon$}
		{
			\Return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]$\;
		}
	}	
%	$D = \emptyset$\;	
}

\caption{Probabilistic $\delta$-reachability \textbf{PHA} and \textbf{NPHA}}
\end{algorithm}

\begin{theorem}\label{thm:dd}
Given $\epsilon \in (0, 1] \cap \mathbb{Q}$, $k\in\mathbb{N}$, and a full {\em\textbf{NPHA}},
there exists an algorithm for computing an interval containing the set of $k$-step reachability 
probabilities. If the system has no nondeterministic parameters, 
the algorithm returns an interval of size not larger than $\epsilon$ containing the
$k$-step reachability probability (\ref{eq:probabilistic-reachability_rv}).
\end{theorem}
Algorithm \ref{alg:gen-prob-reach} drives the whole verification loop, while also handling
discrete random parameters (with essentially the same technique as before). Notice that when
the model has continuous parameters, Algorithm \ref{alg:prob-reach} is utilised.


\begin{algorithm}
\label{alg:gen-prob-reach}
input: continuous random parameters $\overline{r} = \{r_1, \ldots , r_l\}$ with probability 
densities $\overline{f}(r)$, discrete random parameters $\{D_{1}, \ldots, D_p\}$ with probability 
distributions $\overline{p}(\cdot)$, 
$\epsilon_{prod}\in (0,1]\cap\mathbb{Q}$, hybrid system description $\phi$\;
output: interval $[I]$ enclosing reachability probability

\tcp{obtain Cartesian product of discrete random parameters}
$DD$.push$(D_{1}\times \cdots \times D_p)$\;
\tcp{set initial probability intervals for upper/lower approximation}
$P = \emptyset$\;
\While{$size(DD) > 0$}
{
	$\textbf{dd} = DD$.pop()\;
	$m_{\textbf{dd}} = \prod_{d \in \textbf{dd}}p(d)$\;
	\tcp{syntactically replace discrete parameters and get new model}
	\tcp{if no other parameters are present, $\phi_\mathbf{dd}$ becomes an \textbf{HA}}
	$\phi_\textbf{dd} = \phi [(D_1, \ldots D_p) \backslash \textbf{dd}$]\;
	\If{MODEL\_TYPE($\phi_\mathbf{dd}$) = \textbf{HA}}
	{
		\tcp{$\delta$-complete procedure evaluates formula}
		\If{$\phi_\mathbf{dd}$ = $\delta$-sat}
		{
			\tcp{$\delta$-complete procedure evaluates formula}
			\If{$\phi^{C}_\mathbf{dd}$ = $\delta$-sat}
			{
				\tcp{could not correctly decide for $\phi_\mathbf{dd}$}
				$P$.push($m_{\textbf{dd}} \cdot [0.0, 1.0]$)\;
			}\Else
			{
				\tcp{formula $\phi_\mathbf{dd}$ is sat}
				$P$.push($m_{\textbf{dd}} \cdot [1.0, 1.0]$)\;
			}
		} \Else
		{
			\tcp{formula $\phi_\mathbf{dd}$ is unsat}
			$P$.push($m_{\textbf{dd}} \cdot [0.0, 0.0]$)\;
		}
	}
	\If{(MODEL\_TYPE($\phi_\mathbf{dd}$) = \textbf{PHA} $\vee$ \textbf{NPHA})}
	{
		$[I] = \textbf{Algorithm 2}(\overline{r}, \epsilon_{prod}, \phi_{\textbf{dd}})$\;
		\tcp{add obtained probability interval to the stack}
		$P$.push($m_{\textbf{dd}} \cdot [I]$)\;
	}
}
\tcp{obtain sum of all probability intervals on the stack}
$[\underline{[P_{lower}]}, \overline{[P_{upper}]}] = \sum_{[I] \in P}[I]$\;
\Return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]$\;
\caption{Main \preach\ algorithm}
\end{algorithm}


\begin{theorem}\label{thm:complexity}
The complexity of Algorithm \ref{alg:prob-reach} is ${NP}^{(\Sigma_{2}^{P})^{C}}$, where 
$P \subseteq C \subseteq PSPACE$ is the complexity of the terms in the description of the
hybrid system. With Lipschitz-continuous ODEs terms the complexity is $PSPACE$-complete.
\end{theorem}

appendix.tex
\section{Proofs}\label{apndx:proofs}

\begin{proof}[Proposition \ref{prop:BBorel}]
Immediate from the fact that (Definition \ref{def:HS}) the sets defined
by $\text{flow}_q, \text{init}_q, \text{jump}_q$, and $\text{unsafe}_q$ are Borel, and conjunction and
disjunctions correspond to set intersection and union, respectively.
\qed
\end{proof}

\begin{proof}[Proposition \ref{prop:IntegrationComplexity}]
It was proven in \cite[Corollary 6.3]{kobook} that the complexity of computing derivatives of 
$f\in P_{C^5 [a,b]}$ is P. Thus, computing a partial sum on an interval $[x]$ and 
evaluating the formula
\begin{equation*}
width([K]([x])) > \epsilon \frac{width([x])}{b - a}
\end{equation*}
is also polynomial in time. Given an arbitrary partition containing $n$ intervals, the formula above
can thus be verified in polynomial time with respect to the size of the partition. Hence, obtaining a partition such that on each interval the formula above holds is in NP complexity class.
\qed
\end{proof}

We recall that $P_{C^n[a,b]}$ denotes the class of polynomial-time (Type 2) computable functions
whose derivative $f^{(n)}$ exists and is continuous over $[a,b]$. Ko \cite[Section 6.2]{kobook}
showed that if $f$ is also analytic, then integration becomes P. However, such an algorithm
essentially uses truncated Taylor series over an arbitrary partition. Instead, our
Algorithm \ref{alg:integration} adaptively searches for a partition that guarantees the
required error bound $\epsilon$, while having a minimal number of intervals. In practice, this
significantly benefits the performance of our whole implementation. Another advantage of
Algorithm \ref{alg:integration} is that it does not require $f$ to be analytic.


\begin{proof}[Lemma \ref{lemma:converge}]
It was proven in \cite{DBLP:conf/lics/GaoAC12} that satisfiability of a first-order formula implies satisfiability of its weakening. Therefore, following can be equivalently derived:
\begin{equation*}
\phi \rightarrow \phi^{\delta'} \Leftrightarrow \neg \phi \vee \phi^{\delta'} \Leftrightarrow \phi^{\delta'} \vee \neg \phi \Leftrightarrow \neg \phi^{\delta'} \rightarrow \neg \phi
\end{equation*}

Let now $\psi = \phi^{\delta'}$ and $\psi^{\delta^{*}} = \phi^{\delta}$ be weakening of $\psi$. It was proven that if the weakening of the formula is unsatisfiable then the formula is also unsatisfiable. Then:
\begin{equation*}
\neg \psi^{\delta^{*}} \rightarrow \neg \psi \Leftrightarrow \neg \phi^{\delta} \rightarrow \neg \phi^{\delta'}
\end{equation*}
\qed
\end{proof}

\begin{proof}[Proposition \ref{prop:converge}]
By the definition of the decision procedure both formulas can be $\delta$-sat on an interval if and only if the considered interval contains values from the Borel set $B$ and its complement, or when a {\em false alarm} occurs. Then it can be concluded that the initial interval contains a subinterval which is either in the Borel set $B$ or outside it. This can be stated as:
\begin{equation} \label{eq:exists-interval}
\begin{split}
&\exists \delta \in \mathbb{Q}^{+}:(\phi([u, v])-\delta\textbf{-sat}) \wedge (\phi^{C}([u, v])-\delta\textbf{-sat}) \Rightarrow \\ 
&\exists [u', v'] \subseteq [u, v] : ([u', v'] \cap B = [u', v']) \oplus ([u', v'] \cap B = \emptyset)
\end{split}
\end{equation}

Then applying the decision procedure to $[u', v']$ and decreasing $\delta$, it is guaranteed that eventually we will obtain such a $\delta$ that the weakening of the formula will be false. In other words:
\begin{equation} \label{eq:exists-delta-unsat}
\exists \delta \in \mathbb{Q}^{+}: (\neg (\phi([u', v'])-\delta\textbf{-sat})) \oplus (\neg (\phi^{C}([u', v'])-\delta\textbf{-sat}))
\end{equation}

Therefore, by Lemma \ref{lemma:converge} the decision procedure will return \textbf{unsat} for one of the formulas in (\ref{eq:exists-delta-unsat}):
\begin{equation} \label{eq:exists-unsat}
\begin{split}
&\exists \delta \in \mathbb{Q}^{+}: (\neg (\phi([u', v'])-\delta\textbf{-sat})) \oplus (\neg (\phi^{C}([u', v'])-\delta\textbf{-sat})) \Rightarrow \\
&(\phi([u', v'])-\textbf{unsat}) \oplus (\phi^{C}([u', v'])-\textbf{unsat})
\end{split}
\end{equation}
\qed
\end{proof}


\begin{proof}[Proposition \ref{prop:multrv}]

If a hybrid system has $l$ independent initial random parameters with bounded support, then 
the reachability probability can be computed as:
\begin{equation} \label{eq:mulvarind}
  \int_{B} \prod_{i = 1}^{l} dP_i(r_i) = \int_{\Omega} I_{B}(r_{1}, ..., r_{l}) \prod_{i = 1}^{l} dP_i(r_{i})
\end{equation}
where $P_i$ is the probability measure of the $i$-th random parameter $r_i$, $B$ is the 
Borel set (\ref{def:Borel-set}) that contains all the random parameters values for which 
the hybrid system reaches the unsafe region in $k$ steps, $\Omega$ is the domain of the random parameters,
and $I_{B}(r_{1}, ..., r_{l})$ is the indicator function.

In order to compute (\ref{eq:mulvarind}) with precision $\epsilon_{prod}$, we must be able to compute 
\begin{equation} \label{eq:mulvarintegral}
\int_{\Omega} \prod_{i = 1}^{l} dP_i(r_{i})
\end{equation}
with the same precision. By Fubini's theorem, integral (\ref{eq:mulvarintegral}) can be calculated 
as the product 
\[
\int_{\Omega} \prod_{i = 1}^{l} dP_i(r_{i}) = \prod_{i = 1}^{l} \int_{a_{i}}^{b_{i}} dP_i(r_{i}) 
= \prod_{i = 1}^{l} I_{i}
\]
where 
\[
 I_i = \int_{a_{i}}^{b_{i}} dP_i(r_{i}) 
\]
and $a_i, b_i$ are the domain bounds of random parameter $r_i$.

Now, we can compute an interval of length $\epsilon_{i}$ containing the {\em exact} value of each 
integral $I_i$, and let us denote such interval as $[\hat{I}_i, \hat{I}_i +\epsilon_i]$. 
It is thus sufficient to demonstrate how the values $\epsilon_{i}$'s should be chosen 
in order for the integral (\ref{eq:mulvarintegral}) to be contained in an interval of length 
$\epsilon_{prod}$.

According to the rules of interval arithmetics, product of the intervals is contained in the interval:
\begin{equation}\label{eq:intervalprod}
%[\prod_{i = 1}^{l} \hat{I_{i}}, \prod_{i = 1}^{l} \hat{I_{i}} + \epsilon_{prod}] = [\hat{I_{1}}, \hat{I_{1}} 
%+ \epsilon_{1}] \cdot ... \cdot [\hat{I_{l}}, \hat{I_{l}} + \epsilon_{l}] = [\prod_{i = 1}^{l} \hat{I_{i}}, 
%\prod_{i = 1}^{l} (\hat{I_{i}} + \epsilon_{i})]
[\hat{I_{1}}, \hat{I_{1}} + \epsilon_{1}] \cdot [\hat{I_{2}}, \hat{I_{2}} + \epsilon_{2}] \cdots 
[\hat{I_{l}}, \hat{I_{l}} + \epsilon_{l}] 
\ \subseteq\  [\prod_{i = 1}^{l} \hat{I_{i}}, \prod_{i = 1}^{l} (\hat{I_{i}} + \epsilon_{i})]
\end{equation}

Therefore, the $\epsilon_{i}$'s should be chosen such that the interval at the RHS 
of inclusion (\ref{eq:intervalprod}) has length smaller than $\epsilon_{prod}$, \ie, the following should hold:
\begin{equation} \label{eq:mulvarinequality}
\prod_{i = 1}^{l} (\hat{I_{i}} + \epsilon_{i}) - \prod_{i = 1}^{l} \hat{I_{i}} \leq \epsilon_{prod} 
\end{equation}

\iffalse
Let $R$ denote a set $\{I_{1}^{est}, ..., I_{n}^{est}\}$ and $S = \{\epsilon_{1}, ..., \epsilon_{n}\}$. Then the right-hand side of the formula (\ref{eq:mulvarinequality}) can be modified:
\begin{equation}
\begin{split}
\sum ({{R}\choose{0}}{{S}\choose{l}}) + \sum ({{R}\choose{1}}{{S}\choose{l - 1}}) + ... + \sum ({{R}\choose{i}}{{S}\choose{l - i}}) + \\
... + \sum ({{R}\choose{l - 1}}{{S}\choose{1}}) + \sum ({{R}\choose{l}}{{S}\choose{0}})
\end{split}
\end{equation}
where $\sum ({{R}\choose{i}}{{S}\choose{l - i}})$ is a sum of products of all the combinations of length $i$ from the set $R$  multiplied by the combinations of length $l - i$ from set $S$ chosen in such a way that all the indexes in each of the obtained products do not repeat. 

Taking into account that $I_{i}^{est} \in [0, 1]$ the following holds:
\begin{equation} \label{eq:multivarmaininequality}
\begin{split}
\epsilon_{prod} \le \sum {{S}\choose{l}} + \sum {{S}\choose{l - 1}} + ... + \sum {{S}\choose{l - i}} + \\
... + \sum {{S}\choose{1}} + \sum {{S}\choose{0}}
\end{split}
\end{equation}
\fi

Therefore, choosing $\epsilon_{i}$ in such a way that (\ref{eq:mulvarinequality}) holds will guarantee that the {\em exact} value of the product of $l$ integrals is contained in the interval of size $\epsilon_{prod}$. 
If we want all the $\epsilon_{i}$'s equal to a single value $\epsilon$, then formula 
(\ref{eq:mulvarinequality}) can be satisfied by assuming in the worst case $\hat{I}_i = 1$ for all $i$,
which gives
\[
%\epsilon_{prod} \ge {{l}\choose{1}}\epsilon + {{l}\choose{2}}\epsilon^{2} + ... + {{l}\choose{i}}\epsilon^{i} + ... + {{l}\choose{l-1}}\epsilon^{l-1}  + {{l}\choose{l}}\epsilon^{l}
\epsilon_{prod} \geq \prod_{i=1}^l (1+\epsilon) -1 =  \sum_{i=1}^l {{l}\choose{i}} \epsilon^i
\]
where ${{l}\choose{i}}$ is the binomial coefficient.


\qed
\end{proof}

\begin{proposition}\label{prop:bounded}
Given $\epsilon \in (0, 1] \cap \mathbb{Q}$, $k\in\mathbb{N}$ and a hybrid system with one
\textbf{bounded} continuous random initial parameter, there exists an algorithm for computing an interval
of size not larger than $\epsilon$ that contains the value of (\ref{eq:probabilistic-reachability}),
\ie, the probability of reaching the unsafe region in $k$ steps.
\end{proposition}
\begin{proof}%[Proposition \ref{prop:bounded}]
Let $r \in [a, b]$ be a random continuous parameter. Then by applying our validated integration procedure
(Algorithm \ref{alg:integration}) we obtain a partition $\cup_{i = 1}^{n} [r]_{i}$ such that on each of the intervals the value of the partial sum is enclosed by an interval of length $\epsilon \frac{width([r]_{i})}{width([a, b])}$, and the value of the integral on $[a, b]$ is enclosed by the interval of size $\epsilon$.

Let $k$-th step reachability be encoded by the formula $\phi$, and $\phi^{C}$ be derived as in (\ref{eq:phi_C}).
By applying the decision procedure to all the intervals from the initial partition, we can distributed them in three sets $B_{unsat}, B_{C^{unsat}}, B_{\delta-sat}$ containing the intervals where $\phi$ is \textbf{unsat}, $\phi^{C}$ is \textbf{unsat}, and both formulas are $\delta$\textbf{-sat}, respectively. Then the following will hold:
\begin{equation} \label{eq:a-b-integral}
\int_{a}^{b}\, f(r)\, dr = \int_{B_{unsat}}\, f(r)\, dr + \int_{B_{C^{unsat}}}\, f(r)\, dr + \int_{B_{\delta-sat}}\, f(r)\, dr
\end{equation}
The lower and the upper bounds of the interval containing the exact value of the probability can be found as:
\[
\begin{split}
&P_{lower} = \int_{B_{C^{unsat}}}\, f(r)\, dr \\
&P_{upper} = \int_{a}^{b}\, f(r)\, dr - \int_{B_{unsat}}\, f(r)\, dr
\end{split}
\]
Then size of the interval $[P_{lower}, P_{upper}]$ can be calculated as:
\[
\begin{split}
&P_{upper} - P_{lower} = \int_{a}^{b}\, f(r)\, dr - \int_{B_{unsat}}\, f(r)\, dr - \int_{B_{C^{unsat}}}\, f(r)\, dr = \\
&= \int_{B_{\delta-sat}}\, f(r)\, dr
\end{split}
\]

By Proposition \ref{prop:converge} it follows that on each interval $[u, v]$ in $B_{\delta-sat}$ we can obtain a subinterval $[u', v']$ such that it can be added to $B_{unsat}$ or $B_{C^{unsat}}$ and, thus, removed from $B_{\delta-sat}$. Therefore, as $\delta \rightarrow 0$ and $n \rightarrow \infty$ (where $n$ is the number of disjoint subintervals partitioning $B_{\delta-sat}$) the size of set $B_{\delta-sat}$ will be decreasing. Hence, we can conclude that $\int_{B_{\delta-sat}}\, f(r)\, dr \rightarrow 0$, which implies:
\[
\exists \epsilon \in \mathbb{Q}^{+}: P_{upper} - P_{lower} \le \epsilon
\]
\qed
\end{proof}



\begin{proposition} \label{prop:unbounded}
Given $\epsilon \in (0, 1] \cap \mathbb{Q}$, $k\in\mathbb{N}$ and a hybrid system with one
\textbf{unbounded} continuous random initial parameter, there exists an algorithm for computing
an interval of size not larger than $\epsilon$ that contains the value
of (\ref{eq:probabilistic-reachability}).
\end{proposition}
\begin{proof}
Let us recall that calculating the probability of reaching the unsafe region requires integrating an
indicator function with respect the probability measure associated to the random parameter
\begin{equation} \label{eq:probability}
\int_{\Omega}\, I_{B}(r)\, dP(r)
\end{equation}
where $I_{B}(r)$ is the indicator function over set $B$ (\ref{def:Borel-set}), $P$ is the probability 
measure of the random variable, and $\Omega = (-\infty, +\infty)$. In the following we shall simplify
notation and write $dP$ instead of $dP(r)$, since there is only one random variable.

The next inequality can be readily derived from the definition of indicator function:
\[
0 \le \int_{\Omega}\, I_{B}(r)\, dP \le \int_{\Omega}\, dP .
\]
By the property of definite integral, for any $a$ and $b \ge a$:
\[
\int_{\Omega}\, I_{B}(r)\, dP = \int_{a}^{b}\, I_{B}(r)\, dP + \int_{-\infty}^{a}\, I_{B}(r)\, dP + \int_{b}^{\infty}\, I_{B}(r)\, dP .
\]
As $\int_{-\infty}^{a}\, I_{B}(r)\, dP \ge 0$ and $\int_{b}^{\infty}\, I_{B}(r)\, dP \ge 0$, the following 
holds for all $r \in \Omega$
\[
\int_{a}^{b}\, I_{B}(r)\, dP \le \int_{\Omega}\, I_{B}(r)\, dP \le \int_{a}^{b}\, I_{B}(r)\, dP + 1 - \int_{a}^{b}\, dP .
\]
Therefore, the exact value of probability is enclosed by the interval:
\[
\int_{\Omega}\, I_{B}(r)\, dP  \in [\int_{a}^{b}\, I_{B}(r)\, dP, \int_{a}^{b}\, I_{B}(r)\, dP + 1 - \int_{a}^{b}\, dP]
\]

By Proposition \ref{prop:bounded} we can calculate the lower and the upper bounds of the probability 
over the bounded interval $[a, b]$:
\begin{equation} \label{eq:exact-probability}
\int_{\Omega}\, I_{B}(r)\, dP  \in [(\int_{a}^{b}\, I_{B}(r)\, dP)_{lower}, (\int_{a}^{b}\, I_{B}(r)\, dP)_{upper} + 1 - \int_{a}^{b}\, dP]
\end{equation}

Now it is desired that the interval in formula (\ref{eq:exact-probability}) is of length $\epsilon$. For this the error $\epsilon$ can be presented as a sum of two components $\epsilon_{inf}$ and $\epsilon_{prob}$ that are chosen such that: $\epsilon \ge \epsilon_{inf} + \epsilon_{prob}$ where $\epsilon_{inf} \ge 1 - \int_{a}^{b}\, dP$ and $\epsilon_{prob} \ge (\int_{a}^{b}\, I_{B}(r)\, dP)_{upper} - (\int_{a}^{b}\, I_{B}(r)\, dP)_{lower}$.

The values $a$ and $b$ can be obtained by solving the first inequality as a first order formula:
\begin{equation} \label{eq:epsilon-inf}
\exists a \in [u_{a}, v_{a}], \exists b \in [u_{b}, v_{b}] : (\frac{dF}{dx} = f(x)) \wedge (F(a) = 0) \wedge (F(b) \ge 1 - \epsilon_{inf})
\end{equation}
where $f$ is the probability density function of the random parameter, which is known to the user. (Note that
$F$ thus denotes the cumulative distribution function of the random parameter.)
Then the values $a$ and $b$ derived from formula (\ref{eq:epsilon-inf}) are used to compute the interval $[(\int_{a}^{b}\, I_{B}(r)\, dP)_{lower}, (\int_{a}^{b}\, I_{B}(r)\, dP)_{upper}]$ of length $\epsilon_{prob}$. This can be performed for an arbitrary positive rational number (by Proposition \ref{prop:bounded}).

If formula (\ref{eq:epsilon-inf}) is unsatisfiable then it means that bounds for the variables $a$ and $b$ should be enlarged and the formula should be verified again. This process should repeat until the formula is satisfiable and the values $a$ and $b$ are obtained.

\qed
\end{proof}



\begin{proof}[Proposition \ref{prop:cn_nd}]
Evaluating formulas $\phi$ and $\phi^C$ on two boxes $[\overline{r}]$ and $[\overline{z}]$ (over random and nondeterministic continuous parameters, respectively) there are four possible outcomes:
\begin{itemize}
        \item{$\phi([\overline{r}], [\overline{z}])$ is \textbf{unsat}. Hence, there are {\em for sure} no values in $[\overline{r}]$ and $[\overline{z}]$ such that the system reaches the unsafe region, so $[\overline{r}]$ is not in $B$.}
        \item{$\phi([\overline{r}], [\overline{z}])$ is \textbf{$\delta$-sat}. Then, there is a value in $[\overline{r}], [\overline{z}]$ such that the system reaches $U$ or $U^{\delta}$ ($\delta$-weakening of set $U$).}
        \item{$\phi^C([\overline{r}], [\overline{z}])$ is \textbf{unsat}. Therefore, there is {\em for sure} no value in $[\overline{r}]$ and $[\overline{z}]$ such that for all time points on the $k$-th step the system stays within the complement of the unsafe region. In other words, for all the values in $[\textbf{r}]$ the system    reaches the unsafe region, so $[\textbf{r}]$ is fully contained in $B$.}
        \item{$\phi^C([\overline{r}], [\overline{z}])$ is \textbf{$\delta$-sat}. Then there is a value in $[\overline{r}], [\overline{z}]$ such that the system stays within $U^{C}$ or $U^{C^{\delta}}$.}
\end{itemize}
Similarly to the approach used in the proof of Proposition \ref{prop:bounded}, lower and upper bounds of 
the reachability probability can be calculated as:
\[
\begin{split}
&P_{lower} = \int_{B_{C^{unsat}}}\, f(r)\, dr \\
&P_{upper} = \int_{a}^{b}\, f(r)\, dr - \int_{B_{unsat}}\, f(r)\, dr
\end{split}
\]
where $B_{unsat}, B_{C^{unsat}}, B_{\delta-sat}$ containing the boxes where $\phi$ is \textbf{unsat}, $\phi^{C}$ is \textbf{unsat}, and both formulas are $\delta$\textbf{-sat} respectively.
Hence, by refining boxes from $B_{\delta-sat}$ until $max(|[\textbf{r}]|) \le \epsilon$, 
we obtain an interval $[P_{lower}, P_{upper}]$ containing the range of probabilities of reaching 
the unsafe region.
\qed
\end{proof}


\begin{proof}[Theorem \ref{thm:dd}]
Let $\phi$ be a formula describing a hybrid system with discrete random parameters 
$\{{D}_{1}, \ldots, {D}_{p}\}$, and let $p(\cdot)$ denotes their probabilities, \ie,
for each $D_{i}$ we have that $ \sum_{j=1}^{\#D_{i}} p(d_{ij}) = 1$.

Let $\textbf{DD} = {D}_{1} \times \cdots \times {D}_{p}$ be the Cartesian product of discrete parameters.
For each $\textbf{dd} = \{d_{11}, d_{22}, ..., d_{pk}\} \in \textbf{DD}$, 
let $m_{\textbf{dd}} = p(d_{11}) \cdot p(d_{21}) \cdot \cdots \cdot p(d_{pk})$.
Substituting all discrete random parameters with their values from \textbf{dd} we will obtain a hybrid system which can be described by a corresponding formula $\phi_{\textbf{dd}}$.

Now depending of the type of the considered hybrid system we can use one of the algorithms already presented.
\begin{itemize}
        \item{\textbf{PHA} or \textbf{NPHA}: we can apply Algorithm \ref{alg:prob-reach} and obtain a probability interval $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]_{\textbf{dd}}$.}
        \item{\textbf{HA}: we can just use the decision procedure described above and evaluate $\phi_{\textbf{dd}}$ and $\phi^{C}_{\textbf{dd}}$. Then returned value depends on the evaluation outcome:
                \begin{itemize}
                        \item{$\phi_{\textbf{dd}}$-\textbf{unsat} return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]_{\textbf{dd}} = [0.0, 0.0]$}
                        \item{$\phi^{C}_{\textbf{dd}}$-\textbf{unsat} return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]_{\textbf{dd}} = [1.0, 1.0]$}
                        \item{$\phi_{\textbf{dd}}$-\textbf{$\delta-sat$} and $\phi^{C}_{\textbf{dd}}$-\textbf{$\delta-sat$} return $[\underline{[P_{lower}]}, \overline{[P_{upper}]}]_{\textbf{dd}} = [0.0, 1.0]$}
                \end{itemize}}
\end{itemize}

Doing so for each $\textbf{dd} \in \textbf{DD}$ we can obtain the resulting probability interval
\begin{equation*}
[P] = \sum_{\textbf{dd} \in \textbf{DD}}(m_{\textbf{dd}} \cdot [\underline{[P_{lower}]}, \overline{[P_{upper}]}]_{\textbf{dd}})
\end{equation*}
\qed
\end{proof}

\begin{proof}[Theorem \ref{thm:complexity}]
The considered algorithm can be presented as two independent components: validated integration and probability calculation.

The decision procedure used in the algorithm consists of two formulas: $\phi$ and $\phi^{C}$,
which are $\Sigma_{1}$ and $\Sigma_{2}$ sentences. Solving these formulas as a $\delta$-SMT problem
is in $(\Sigma_{1}^{P})^{C}$ and $(\Sigma_{2}^{P})^{C}$ complexity classes respectively,
where $P \subseteq C \subseteq PSPACE$ is the complexity of the terms in the 
formula \cite{DBLP:conf/lics/GaoAC12}. 
Hence, the decision procedure is in $(\Sigma_{2}^{P})^{C}$.
Then verification of an arbitrary partition of $n$ intervals is also in $(\Sigma_{2}^{P})^{C}$.
Hence, it is clear that obtaining the {\em correct} partition (such that 
$\overline{[P_{upper}]} - \underline{[P_{lower}]} \le \epsilon_{prob}$) is in 
$NP^{(\Sigma_{2}^{P})^{C}}$.
By Proposition \ref{prop:IntegrationComplexity} the complexity of the verified integration is $NP$,
which is in $NP^{(\Sigma_{2}^{P})^{C}}$. The complexity of the whole algorithm is thus 
$NP^{(\Sigma_{2}^{P})^{C}}$, where $P \subseteq C \subseteq PSPACE$.

Finally, it has been shown in \cite{DBLP:conf/lics/GaoAC12} that if $\phi$ (and thus $\phi^C$) 
includes Lipschitz-continuous ODEs then the $\delta$-SMT problem becomes $PSPACE$-complete.
This lifts the complexity of the whole algorithm to $PSPACE$-complete.

\qed
\end{proof}

\section{Models}\label{apndx:Models}
We give here more information about the models used for our experiments.

\subsection{2D-moving Bouncing ball}
The ball is launched from position $(S_{x} \in [-5, 5], S_{y} = 0)$ with initial speed 
$\upsilon_0 \sim N(20,1)$, \ie, normal distribution with mean 20 and variance 1, and angle 
$\alpha$ to horizon (measured in radians) with the following probability distribution: 
\[
\begin{array}{l}
P[\alpha = 0.5236] = 0.3\\
P[\alpha = 0.7854] = 0.5\\
P[\alpha = 1.0472] = 0.2
\end{array}
\]
After each jump the speed of the ball is multiplied by 0.9. The gravity of Earth parameter
$g\in[9.8,9.81]$ is also nondeterministic. The system is modelled as a hybrid system with 
one mode with dynamics governed by a system of ODEs: 
\[
\begin{array}{l}
S'_{x}(t) = \upsilon_0 \cos{\alpha} \\[1ex]
S'_{y}(t) = \upsilon_0 \sin{\alpha} - gt
\end{array}
\]
The goal of the experiment is to calculate the probability of reaching the region 
$S_{x}(t) \ge 100$ within 0 and 1 jump. The results are presented in Table \ref{table:bouncing-ball}. 
Monte Carlo simulation of continuous nondeterminism in MATLAB was achieved as explained
in Section \ref{sec:Exp}, using uniform discretisation of the domains of the nondeterministic 
parameters ($S_{x}$ and $g$ were discretised with 100 and 10 values, respectively). 
In Figure \ref{fig:bb-graph0} and \ref{fig:bb-graph1} we plot the Monte Carlo reachability 
probability estimate with respect to the nondeterministic parameters $S_x(0)$ and $g$, for 0 and 1 jump,
respectively.
\begin{figure}[ht!] 
	\centering
	\includegraphics[width=120mm]{bb-k-0.png}
	\caption{Monte Carlo simulation of the bouncing ball model: Reachability probability ($k=0$) 
	estimate $P$ with respect to nondeterministic parameters $S_{x}(0)\in [-5,5]$ and $g\in [9.8,9.81]$.}
	\label{fig:bb-graph0}
\end{figure}

\begin{figure}[ht!] 
	\centering
	\includegraphics[width=120mm]{bb-k-1.png}
	\caption{Monte Carlo simulation of the bouncing ball model: Reachability probability ($k=1$) 
	estimate $P$ with respect to nondeterministic parameters $S_{x}(0)\in [-5,5]$ and $g\in [9.8,9.81]$.}
	\label{fig:bb-graph1}
\end{figure}


\begin{figure}[ht!] 
	\centering
	\includegraphics[width=90mm]{bouncing-ball-k-1-random-speed.png}
	\caption{2D-moving bouncing ball scenario}
	\label{fig:bouncing-ball}
\end{figure}

\begin{table}[ht!]
\begin{adjustwidth}{-1in}{-1in}
\centering
\caption{2D-moving bouncing ball model; see legend in Section \ref{sec:Exp}.}
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method & Model type& $k$ & $\epsilon$ & $length$ & Probability interval & $CPU$ & \\ [0.5ex]
\hline
\multirow{2}{4em}{{\sffamily Prob Reach}} &\textbf{NPHA}& 0 & $10^{-3}$ & $2.38\cdot10^{-4}$ & [0.000013103, 0.000250681] & 223 & \\ [0.5ex]
&\textbf{NPHA}& 1 & $10^{-3}$ & $6.464\cdot10^{-2}$ &[0.0647381, 0.12937951] & 1,605 & \\ [0.5ex]
\hline \hline
Method & Model type & $k$ & $\zeta$ & $c$ & Monte Carlo interval & $CPU$ & $N$\\ [0.5ex]
\hline
\multirow{2}{4em}{Monte Carlo}&\textbf{NPHA}& 0 & $5\cdot 10^{-3}$ & 0.99 & [0, 0.00520629] & 1,482 & 92,104\\ [0.5ex]
&\textbf{NPHA}& 1 & $5\cdot10^{-3}$ & 0.99 & [0.0585, 0.1367] & 1,485 & 92,104 \\ [0.5ex]
\hline
\end{tabular}
\label{table:bouncing-ball}
\end{adjustwidth}
\end{table}


\hide{
\begin{table}[ht] 
\caption{Computing probabilistic reachability for the bouncing ball model}
\centering
\begin{tabular}{c c c c r}
\hline\hline \\ [0.5ex]
$k$ & $\delta$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
0 & $10^{-9}$ & [8.2175709429082528e-05, 8.2176188377012998e-05] & 64 & 7\\ [0.5ex]
1 & $10^{-9}$  & [0.13794836312721623, 0.13794836411185463] & 192 & 29\\ [0.5ex] 
2 & $10^{-9}$  & [0.50868960502622063, 0.50868960601267332] & 927 & 164\\ [0.5ex] 
3 & $10^{-9}$  & [0.73876740050846257, 0.73876740135236485] & 3806 & 563\\ [0.5ex] 
\hline \\ [0.5ex]
\end{tabular} 

$k$ = number of discrete transitions, $\delta$ = size of probability interval, $CPU_{seq}$ = CPU time of sequential version of the algorithm in seconds, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
\label{table:bouncing-ball} 
\end{table} 

\begin{table}[ht!] 
\caption{Monte Carlo simulations for the bouncing ball model}
\centering
\begin{tabular}{c c c r}
\hline\hline \\ [0.5ex]
$k$ & $P$ & $CI$ & $CPU$\\ [0.5ex] 

\hline \\ [0.5ex]
0 & 8.214014e-05 & [7.714014e-05, 8.714014e-05] & 14,394\\ [0.5ex]
1 & 0.137946 & [0.137941, 0.137951] & 16,531\\ [0.5ex] 
2 & 0.5086959 & [0.5086909, 0.5086959] & 18,052\\ [0.5ex] 
3 & 0.7387694 & [0.7387644, 0.7387744] & 17,530\\ [0.5ex] 
\hline \\ [0.5ex]
\end{tabular} 

$k$ = number of discrete transitions, $P$ = probability estimate, $CI$ = confidence interval. Sample size = 230,258,509,300.
\label{table:bouncing-ball-MATLAB} 
\end{table}
}

\iffalse
\subsection{Thermostat} 
The main purpose of the system is to keep the temperature within the desired range. 
The system is modelled by a two mode hybrid system \cite{DBLP:conf/hybrid/AlurCHH92} 
(Fig. \ref{fig:thermostat-2m}). The temperature is changing exponentially and it is 
decreasing in the first mode and increasing in the second mode.
The system starts in mode 1 with the initial temperature $T_{0}$ which is normally 
distributed ($\mu = 30$ and $\sigma = 1$). When the temperature drops to the minimum 
level $T_{min} = 18$, the system makes a transition to mode 2, where the temperature 
increases until it reaches a maximum level $T_{max} = 22$. Then the system makes a 
jump to mode 1 and the loop repeats again. In the model we use function $\tau(t)$ 
to represent the global time, as the current time is reset when the system makes 
a discrete transition.

\begin{figure}[ht!] 
\centering
\includegraphics[width=90mm]{thermostat-2m}
\caption{A 2-mode thermostat hybrid system}
\label{fig:thermostat-2m}
\end{figure}

The goal of the experiment is to calculate the probability of reaching the region 
$T(t) \in [19.9, 20.1]$ in mode 2 at various time points $\tau$. The results are 
presented in Table \ref{table:thermostat-2m}.

\begin{table}[ht] 
\caption{Computing probabilistic reachability of the thermostat model}
\centering
\begin{tabular}{c c c c c c}
\hline\hline \\ [0.5ex]
\# & $k$ & $\tau$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
1 & 1 & 0.6 & [0.0066784445550149884, 0.0066784455059296028] & 71 & 7\\ [0.5ex]
2 & 5 & 1.8 & [0.0026170599124024315, 0.0026170608351848692] & 213 & 23\\ [0.5ex] 
3 & 7 & 2.4 & [0.0015794358437303552, 0.0015794367918444926] & 364 & 49\\ [0.5ex] 
\hline \\ [0.5ex]
\end{tabular} 
\label{table:thermostat-2m}

$k$ = number of discrete transitions, $\tau$ = global time, $CPU_{seq}$ = CPU time of sequential version of the algorithm in seconds, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
{\em CPU} = CPU time in seconds

\end{table}

We can extend the 2-mode thermostat model to a 4-mode version by adding two 
delay modes (Fig. \ref{fig:thermostat-4m}).
\begin{figure}[ht!] 
\centering
\includegraphics[width=90mm]{thermostat-4m}
\caption{A 4-mode thermostat hybrid system}
\label{fig:thermostat-4m}
\end{figure}
The initial mode of the system is mode 1. Modes 1 and 3 are equivalent to modes 1 and 2 in 
the 2-mode thermostat. Modes 2 and 4 model a delay of 0.1 seconds.
The goal of the experiment is to calculate the probability of reaching the region 
$T(t) \in [19.9, 20.1]$ in mode 3 at various time points $\tau$. The results are presented in 
Table \ref{table:thermostat-4m}.

\begin{table}[ht!] 
\caption{Probabilistic reachability of the thermostat model (4 modes)}
\centering
\begin{tabular}{c c c c c c}
\hline\hline \\ [0.5ex]
\# & $k$ & $\tau$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
1 & 2 & 0.6 & [0.0, 8.5404844380251995e-11] & 52 & 4\\ [0.5ex]
2 & 6 & 1.7 & [9.4398656477534736e-8, 9.5194741708276046e-8] & 268 & 28\\ [0.5ex] 
3 & 6 & 1.8 & [0.0039559433434726825, 0.0039559441734411333] & 578 & 75\\ [0.5ex] 
\hline \\ [0.5ex]
\end{tabular} 
\label{table:thermostat-4m} 

$k$ = number of discrete transitions, $\tau$ = global time, $CPU_{seq}$ = CPU time of sequential version of the algorithm in seconds, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
{\em CPU} = CPU time in seconds

\end{table}
\fi

\hide{
\subsection{Controlled bouncing ball}
Consider a 2-mode hybrid system (Fig. \ref{fig:controlled-bouncing-ball}) modelling a controlled 
bouncing ball \cite{ADHS09}. In mode 1, a ball of mass $m = 7$ is dropped on a platform attached 
to a stiff spring and a damper from a random height $H_{0}$, which is distributed normally 
($\mu = 9$ and $\sigma = 1$). When the ball reaches the platform ($H = 0$) the system makes a 
transition to mode 2, where the ball is reflected from the platform and it jumps back to mode 1 
when the height of the ball is greater than 0.

\begin{figure}[ht!] 
\centering
\raisebox{-0.5\height}{\includegraphics[width=40mm]{cbb-picture}}
\hspace*{.2in}
\raisebox{-0.5\height}{\includegraphics[width=70mm]{controlled-bouncing-ball}}
\caption{A figure (left hand side) and a model (right hand side) of a controlled bouncing ball with $R = 5$, $C = 0.0025$ and $g = 9.8$}
\label{fig:controlled-bouncing-ball}
\end{figure}

The goal of the experiment is to calculate the probability that the ball reaches 
the region $H >= 7$ in mode 1 after making one bounce. The results are presented 
in Table \ref{table:controlled-bouncing-ball}.


\iffalse
\begin{table}[ht!] 
\caption{Probabilistic reachability in the controlled bouncing ball model}
\centering
\begin{tabular}{c c c c c}
\hline\hline \\ [0.5ex]
\# & $\epsilon$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
1 & $10^{-2}$ & [0.1999606932602, 0.2071983644358] & 70 & 15\\ [0.5ex]
2 & $10^{-9}$ & [0.2049030217063, 0.2049030224545] & 8332 & 2581\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:controlled-bouncing-ball} 

$\epsilon$ = width of the probability interval, $CPU_{seq}$ = CPU time of sequential version of the algorithm in seconds, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
\end{table}

The obtained result was validated using Monte Carlo simulation in MATLAB 
(Table \ref{table:contolled-bouncing-ball-MATLAB}). The time variable was discretised 
over the interval $[0, 5]$ in mode 1 and $[0, 0.5]$ in mode 2, 
and the system was simulated using random samples distributed normally with mean 
$\mu = 9$ and standard deviation $\sigma = 1$. The number of samples was calculated using 
the Chernoff-Hoeffding bound with $\epsilon = 0.005$ and $c = 0.99$,
which returned $N = 92,104$.

\begin{table}[ht!] 
\caption{Monte Carlo simulation for the controlled bouncing ball model}
\centering
\begin{tabular}{c c c}
\hline\hline \\ [0.5ex]
$P$ & $CI$ & $CPU$\\ [0.5ex] 

\hline \\ [0.5ex]
0.2045948 & [0.1995948, 0.2095948] & 50,528\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 

$P$ = probability estimate, $CI$ = confidence interval and {\em CPU} = CPU time in seconds. Sample size = 92,104.
\label{table:contolled-bouncing-ball-MATLAB} 
\end{table}
\fi
}

\subsection{Starvation model}
In humans, enduring fasting for 3-4 days will consume all the glucose reserves of the body.
At this point, the energy to sustain the human body is produced from fat $F(t)$, muscles $M(t)$ 
and ketone bodies $K(t)$ (for brain function) \cite{Starvation}. The ODE system below represents 
the dynamics of the described variables:
\begin{equation*}
\begin{split}
\frac{dF}{dt} &= F(\frac{- a}{1 + K} - \frac{1}{\lambda_{F}}  (\frac{C + g L_{0}}{F + M} + g)) \\
\frac{dM}{dt} &= - \frac{M}{\lambda_{M}}  (\frac{C + g L_{0}}{F + M} + g) \\
\frac{dK}{dt} &= \frac{V a F}{1 + K} - b
\end{split}
\end{equation*}
We consider two scenarios where parameter $g\sim N(10.96,1)$, \ie, normally distributed with
mean 10.96 and variance 1, and:
\begin{itemize}
	\item{$b \in [0.05, 0.075]$ is nondeterministic}; or
	\item{$b$ is a discrete random parameter with the probability distribution: }
\[
\begin{array}{l}
P[b = 0.05] = 0.1\\
P[b = 0.06] = 0.2\\
P[b = 0.07] = 0.3\\
P[b = 0.075] = 0.4
\end{array}
\]
\end{itemize}

The probabilistic reachability property investigated in the experiment is: {\em what is 
the probability that muscle mass will decrease by 40\% within 25 days?}
Numerical values for all deterministic parameters in the model are presented 
in Table \ref{table:starvation-param} and verification results are featured in 
Table \ref{table:starvation}. Monte Carlo simulation of continuous nondeterminism 
in MATLAB was achieved as explained in Section \ref{sec:Exp}, via uniform discretisation 
(10 values) of the nondeterministic parameter $b$. In Figure \ref{fig:starvation-graph} we
plot the Monte Carlo reachability probability and confidence interval with respect to the 
value of parameter $b$.

\begin{table}[ht!] 
\caption{Starvation model parameters and initial conditions}
\centering
\begin{tabular}{c c c c c c c c c c}
\hline\hline 
Param. & Value & Param. & Value & Param. & Value & Param. & Value & Param. & Value\\ [0.5ex] 

\hline 
$a$ & 0.013 & $\lambda_{F}$ & 7777.8 & $M(0)$ & 43.6 & $V$ & 0.9 & $K(0)$ & 0.02\\ [0.5ex]
$C$ & 772.3 & $\lambda_{M}$ & 1400 & $F(0)$ & 25 & $L_{0}$ & 30.4 &  & \\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:starvation-param} 

\end{table}

\begin{figure}[ht!] 
\centering
\includegraphics[width=120mm]{starvation-graph.png}
\caption{Monte Carlo simulation of the starvation model: Reachability probability estimate $P$ 
(solid line) with respect to nondeterministic parameter $b\in [0.05, 0.075]$. For each 
(discretised) value of $b$  we give a Chernoff-Hoeffding confidence interval, denoted by dotted lines.}
\label{fig:starvation-graph}
\end{figure}

\iffalse
\begin{table}[ht!] 
\caption{Probabilistic reachability in the starvation model}
\centering
\begin{tabular}{c c c c}
\hline\hline \\ [0.5ex]
$\epsilon$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
$10^{-2}$ & [0.919528, 0.927383] & 160 & 5\\ [0.5ex]
$10^{-3}$ & [0.921598, 0.922369] & 53 & 8\\ [0.5ex]
$10^{-9}$ & [0.9221401013279, 0.9221401020929] & 298 & 41\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:starvation} 

$\epsilon$ = width of the probability interval, $CPU_{seq}$ = CPU time of sequential version of the algorithm in seconds, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
\end{table}

The obtained result was validated using Monte Carlo simulation in MATLAB 
(Table \ref{table:starvation}). The time variable was discretised 
over the interval $[0, 25]$, and the system was simulated using random samples distributed normally with mean 
$\mu = 10.96$ and standard deviation $\sigma = 1$. The number of samples was calculated using 
the Chernoff-Hoeffding bound with $\epsilon = 0.005$ and $c = 0.99$,
which returned $N = 92,104$.

\begin{table}[ht] 
\caption{Monte Carlo simulation for the starvation model}
\centering
\begin{tabular}{c c c}
\hline\hline \\ [0.5ex]
$P$ & $CI$ & $CPU$\\ [0.5ex] 

\hline \\ [0.5ex]
0.9211543 & [0.9161543, 0.9261543] & 36,396\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 

$P$ = probability estimate, $CI$ = confidence interval and {\em CPU} = CPU time in seconds. Sample size = 92,104.
\label{table:starvation-MATLAB} 
\end{table}
\fi


\subsection{Road scenario}

\begin{figure}[ht!] 
\centering
\includegraphics[width=120mm]{cars.png}
\caption{Road scenario}
\label{fig:cars}
\end{figure}

We consider the road scenario inspired by a model presented in \cite{SMCnd} and depicted in 
Figure \ref{fig:cars}. Two cars ($Car\#1$ and $Car\#2$) move
on the same lane, starting at coordinates $S_{01} = 0$ and $S_{02} = S_{01} + \upsilon_{1} \cdot t_{safe}$,
where $t_{safe} = 2 sec$ implements the so-called ``two seconds rule'' for maintaining a safety 
distance between two cars.

We describe a car collision scenario and we model it with the hybrid automaton given in 
Figure \ref{fig:cars_model}. Starting in Mode 1 at time $t = 0$, $Car\#1$ changes lane and starts 
accelerating at $a_{a1}$ $m/s^{2}$, while $Car\#2$ is moving in the initial lane 
with speed $\upsilon_{2}$. Upon reaching the maximum speed $\upsilon_{max}$, the
system switches to Mode 2, where $Car\#1$ keeps moving at this speed until it gets ahead 
of $Car\#2$ by the safety distance $S_{safe} = \upsilon_{2} \cdot t_{safe}$.
After that, we switch to Mode 3: $Car\#1$ returns to the initial lane and starts decelerating 
at $a_{d1}$. For the driver of $Car\#2$ it takes $t_{react} = 1 sec$ to react
(the system switches to Mode 4) and then it starts decelerating as well (with random acceleration 
$a_{d2} \sim N(-1.35, 0.01)$). In Mode 3, 4, and 5 we also have an invariant specifying that $Car\#1$
should precede $Car\#2$ at all time. We calculate the {\em probability of observing a car collision 
in Mode 5}, where $Car\#1$ is stopped. Numerical values for all deterministic parameters in the model 
are given in Table \ref{table:cars-param} are verification results are presented in Table \ref{table:cars}.

\begin{table}[ht!] 
\caption{Car collision model parameters and initial conditions}
\centering
\begin{tabular}{c c c c c c c c c c}
\hline\hline 
Param. & Value & Param. & Value & Param. & Value & Param. & Value & Param. & Value\\ [0.5ex] 

\hline 
$\upsilon_{1}$ & 11.12 & $\upsilon_{2}$ & 11.12 & $\upsilon_{max}$ & 16.67 & $t_{safe}$ & 2 & $S_{safe}$ & $\upsilon_{2} \cdot t_{safe}$\\ [0.5ex]
$a_{a1}$ & 3 & $a_{d1}$ & -4 & $t_{react}$ & 1 & $S_{01}$ & 0 & $S_{02}$ & $S_{01} + \upsilon_{1} \cdot t_{safe}$ \\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:cars-param} 

\end{table}

\begin{figure}[ht!] 
\centering
\includegraphics[width=120mm]{cars_model.pdf}
\caption{Car collision model}
\label{fig:cars_model}
\end{figure}


\subsection{Prostate cancer therapy}
We consider a model of personalised prostate cancer therapy introduced by 
Ideta {\em et al.}~\cite{Ideta2008} and improved by Liu {\em et al.}~\cite{Liu2015}. 
Intermittent androgen suppression (IAS) has proved to be more effective than constant androgen 
suppression (CAS) in delaying the recurrence of prostate cancer. Briefly, the personalised 
therapy comprises of two repeating stages. The patient's prostate-specific antigen (PSA) level 
is monitored throughout the therapy. When the PSA level reaches an upper threshold, the patient 
starts receiving treatment ({\em on-therapy} stage) until the PSA level decreases to 
a lower threshold ({\em off-therapy}). The main aim of the therapy is to delay cancer relapse 
for as long as possible.

The model of the therapy is given in Figure \ref{fig:cancer_model} (a full explanation of the 
model and its parameters can be found in \cite{Liu2015}). Mode 1 is the {\em on-therapy} 
stage, and it continues until the PSA level (measured by $x+y$) is above threshold $r_{0} = 4$. 
Then the system makes a transition to the {\em off-therapy} mode which continues until PSA level 
is below $r_{1} = 10$. We explore the following scenarios:
\begin{itemize}
	\item{$\alpha_{y}$ is distributed normally ($N(0.05,0.01)$)} and $\alpha_{x}= 0.0197$
	\item{$\alpha_{y}$ is distributed normally ($N(0.05,0.01)$) and $\alpha_{x} \in [0.0197,0.0204]$ is nondeterministic}
\end{itemize}
For the cases above we calculate the probability of cancer relapse (\ie, $y\geq 1$) within 100 
days of using the personalised cancer therapy. Numerical values of all the parameters in the model 
are presented in the Table \ref{table:cancer-param} and verification results are featured in 
Table \ref{table:cancer}. Monte Carlo simulation of continuous nondeterminism in MATLAB 
was achieved as detailed in Section \ref{sec:Exp}, using uniform discretisation (20 values) of
the domain of the nondeterministic parameter $\alpha_{x}$. In Figure \ref{fig:cancer-graph} we
plot Monte Carlo reachability probability and confidence interval with respect to the value 
of parameter $\alpha_x$.

\begin{figure}[ht!] 
\centering
\includegraphics[width=120mm]{cancer-graph.png}
\caption{Monte Carlo simulation of the prostate cancer therapy model: Reachability probability estimate
$P$ (solid line) with respect to the nondeterministic parameter $\alpha_x\in [0.0197,0.0204]$. 
For each (discretised) value of $\alpha_x$ we give a Chernoff-Hoeffding confidence interval, 
denoted by dotted lines.}
\label{fig:cancer-graph}
\end{figure}

\begin{table}[ht!] 
\caption{Prostate cancer therapy model parameters and initial conditions}
\centering
\begin{tabular}{c c c c c c c c c c}
\hline\hline 
Param. & Value & Param. & Value & Param. & Value & Param. & Value & Param. & Value\\ [0.5ex] 

\hline 
$\beta_{x}$ & 0.0175 & $\beta_{y}$ & 0.0168 & $k_1$ & 10.0 & $k_2$ & 1.0 & $k_3$ & 10.0\\ [0.5ex]
$k_4$ & 2 & $m_1$ & $10^{-5}$ & $z_0$ & 12 & $\gamma$ & 0.08 & $r_1$ & 10.0\\ [0.5ex]
$r_0$ & 4.0 & $d_0$ & $1.0$ & $c_1$ & 0.01 & $c_2$ & 0.03 & $c_3$ & 0.02\\ [0.5ex]
$x(0)$ & 19 & $y(0)$ & $0.1$ & $z(0)$ & 12.5 &  &  &  & \\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:cancer-param} 
\end{table}

\begin{figure}[] 
\centering
\includegraphics[width=120mm]{cancer_model.pdf}
\caption{Personalized prostate cancer therapy model}
\label{fig:cancer_model}
\end{figure}


\subsection{Prostate cancer therapy: \preach\ file}
\begingroup
\fontsize{7pt}{7pt}\selectfont
\begin{Verbatim}[frame=single, numbers=left, samepage = true]
// This is a pdrh file corresponding to the prostate cancer therapy model with 
// one random and one nondeterministic parameter.
MODEL_TYPE(NPHA) // defining model type
#define betax 0.0175
#define betay 0.0168
#define k1 10.0
#define k2 1.0
#define k3 10.0
#define k4 2
#define m1 0.00005
#define z0 12.0
#define gamma 0.08
#define r1 10.0
#define r0 4.0
#define d0 1.0
#define c1 0.01
#define c2 0.03
#define c3 0.02
#define Gx ((alphax/(1+exp((k1-z)*k2)))-(betax/(1+exp((z-k3)*k4))))
#define Gy ((alphay * (1 - (d0 * (z / z0)))) - betay)
#define Mxy (m1 * (1 - (z / z0)))
#define scale 1.0
#define T 100.0
N(0.05,0.01)alphay; // random parameter, normally distributed
[0,T]time;
[0,T]tau;
[0,100.0]x;
[0,10.0]y;
[0.0,100.0]z;
[0.0197,0.0204]alphax; // nondeterministic parameter
{ 
mode1; // on-therapy
	invt:
		(y <= 1);
	flow:
		d/dt[x]=scale * ((Gx - Mxy - c1) * x + c2);
		d/dt[y]=scale * (Mxy * x + Gy * y);
		d/dt[z]=scale * (-z * gamma + c3);
		d/dt[tau]=scale * 1.0;
	jump:
		((x+y)=r0)==>@2(and(tau'=tau)(x'=x)(y'=y)(z'=z));
}
{ 
mode2; // off-therapy
	invt:
		(y <= 1);
	flow:
		d/dt[x]=scale * ((Gx - Mxy - c1) * x + c2);
		d/dt[y]=scale * (Mxy * x + Gy * y);
		d/dt[z]=scale * ((z0 - z) * gamma + c3);
		d/dt[tau]=scale * 1.0;
	jump:
		((x+y)=r1)==>@1(and(tau'=tau)(x'=x)(y'=y)(z'=z));
}
init:
	@1(and (x = 19) (y = 0.1) (z = 12.5) (tau = 0));
goal: // unsafe region
	@2(and(y <= 1)(tau = T));
goal_c: // unsafe region complement
	@2(and(y > 1.0)(tau < T));
\end{Verbatim}
\endgroup




\hide{
\subsection{Insulin-glucose regulatory model}
We consider an insulin-glucose regulatory system for patients with type-1 diabetes. 
Our tool was applied to the insulin infusion model introduced in 
\cite{Sankaranarayanan:2012:SII:2415548.2415569} based on Hovorka {\em et al.}'s glucoregulatory 
model \cite{Hovorka02,Hovorka04}. The model consists of four subsystems: meal, insulin pump, 
glucoregulatory model and monitor (Figure \ref{fig:insulin-infusion-high-level}). 

\begin{figure}[ht!] 
\centering
\includegraphics[width=120mm]{insulin-infusion-high-level.pdf}
\caption{Insulin-glucose regulatory high level model}
\label{fig:insulin-infusion-high-level}
\end{figure}

The meal parameters (\eg, amount of carbohydrates in the the meal and their glycemic index) 
are introduced into the meal subsystem by the patient prior to meal consumption. It models 
the food absorption by the human guts and outputs a glucose absorption rate. According to the meal
characteristics of the patient the system calculates and inputs the initial insulin infusion rate.
These two parameters are used by the glucoregulatory subsystem evaluating a glucose level of the
patient which is used as an input parameter for the monitor. When the glucose level goes outside
the predefined corridor the monitor sends a signal to the pump to increase or decrease the insulin
infusion rate. We consider a simplified version of the glucoregulatory system assuming that glucose 
absorption rate ($UG$) and insulin infusion rate ($u_0$) are given explicitly. The ODEs 
are presented in (\ref{eq:insulin-infusion-low-level}); the model parameters are given
in Table \ref{table:insulin-glucose-param}.

\begin{equation}\label{eq:insulin-infusion-low-level}
\begin{split} 
\frac{dQ_{1}}{dt} &= -F^{c}_{01} - x_{1}Q_{1} + k_{12}Q_{2}-F_{R}+EGP_{0}(1-x_{3})+ 0.18UG\\
\frac{Q_{2}}{dt} &=x_{1}Q_{1}-(k_{12} + x_{2})Q_{2}\\
\frac{dS_{1}}{dt} &=u - \frac{S_{1}}{t_{maxI}}\\
\frac{dS_{2}}{dt} &=\frac{S_{1}-S_{2}}{t_{maxI}}\\
\frac{dI}{dt} &=\frac{S_{2}}{t_{maxI}V_{I}}-k_{e}I\\
\frac{dx_{1}}{dt} &=-k_{a1}x_{1}+k_{b1}I\\
\frac{dx_{2}}{dt} &=-k_{a2}x_{2}+k_{b2}I\\
\frac{dx_{3}}{dt} &=-k_{a3}x_{3}+k_{b3}I\\
F_{01}^{c} &=\frac{F_{01}G}{0.85(G+1)}\\
G &=\frac{Q_{1}}{V_{G}}
\end{split}
\end{equation}

We consider the following scenario. Initially a meal is consumed ($UG = 8$) and insulin pump does not infuse any insulin ($u_{0}=0$). The glucose level starts rising and when it reaches the point of $G = 10$ the monitor sends a signal to the pump to increase the infusion rate ($u = 0.36$). Randomising $x_{3}(0)$ normally ($\mu = 0.05$ and $\sigma = 0.005$) we want to calculate the probability that the glucose level returns back to normal within 60 minutes after the pump started infusion. The results obtained using our tool are given in Table \ref{table:insulin-glucose}.


\begin{table}[ht!] 
\caption{Insulin-glucose regulatory model parameters and initial conditions}
\centering
\begin{tabular}{c c c c c c c c c c}
\hline\hline 
Param. & Value & Param. & Value & Param. & Value & Param. & Value & Param. & Value\\ [0.5ex] 
\hline 
$Q_{1}(0)$ & 64.0 & $S_{1}(0)$ & 4.2 & $I(0)$ & 0.03 & $x_{2}(0)$ & 0.045 & $V_{I}$ & $0.12w$\\ [0.5ex]
$Q_{2}(0)$ & 40.0 & $S_{2}(0)$ & 4.0 & $x_{1}(0)$ & 0.03 & $w$ & 100 & $V_{G}$ & $0.16w$\\ [0.5ex]
$k_{a1}$ & 0.006 & $k_{a2}$ & 0.06 & $k_{a3}$ & 0.03 & $k_{b1}$ & 0.0034 & $k_{b2}$ & 0.056\\ [0.5ex]
$k_{b3}$ & 0.024 & $k_{e}$ & 0.138 & $k_{12}$ & 0.066 & $t_{max,I}$ & 55 & $F_{01}$ & $0.0097w$\\ [0.5ex]
$F_{R}$ & 0.0 & $EGP_{0}$ & $0.0161w$ & $u_{0}$ & 0.0 & $u$ & 0.36 & $UG$ & 8\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:insulin-glucose-param} 

\end{table}





\iffalse
\begin{table}[ht!] 
\caption{Probabilistic reachability in the insulin-glucose regulatory model}
\centering
\begin{tabular}{c c c c c}
\hline\hline \\ [0.5ex]
$\epsilon$ & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline \\ [0.5ex]
$10^{-2}$ & [0.994589, 0.999917] & 2,805,634 & 165,404\\ [0.5ex]
$10^{-3}$ & [0.999107, 0.999917] & 3,326,581 & 443,910\\ [0.5ex]
$10^{-4}$ & [0.999657, 0.999712] & 3,498,765 & 490,257\\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 
\label{table:insulin-glucose} 

$\epsilon$ = width of the probability interval, $CPU_{par}$ = CPU time of parallel version of the algorithm in seconds
\end{table}
\fi

The obtained result was validated using Monte Carlo simulation in MATLAB 
(Table \ref{table:insulin-glucose}). The time variable was discretised 
over the interval $[0, 200]$, and the system was simulated using random 
samples distributed normally with mean $\mu = 0.05$ and standard deviation 
$\sigma = 0.005$. The number of samples was calculated using 
the Chernoff-Hoeffding bound with $\zeta = 0.0025$ and $c = 0.99$,
which returned $N = 92,104$.

\iffalse
\begin{table}[ht] 
\caption{Monte Carlo simulation for the insulin-glucose regulatory model}
\centering
\begin{tabular}{c c c}
\hline\hline \\ [0.5ex]
$P$ & $CI$ & $CPU$\\ [0.5ex] 

\hline \\ [0.5ex]
0.99956 & [0.99706, 1.00206] & 219,623 \\ [0.5ex]
\hline \\ [0.5ex]
\end{tabular} 

$P$ = probability estimate, $CI$ = confidence interval and {\em CPU} = CPU time in seconds. Sample size = 92,104.
\label{table:insulin-glucose-MATLAB} 
\end{table}
\fi
}


\section{$\delta$-satisfiability}\label{apndx:delta-sat}

In order to overcome the undecidability of reasoning about general real formulae, 
Gao {\em et al.} recently
defined the concept of $\delta$-satisfiability over the reals \cite{DBLP:conf/lics/GaoAC12}, 
and presented a corresponding $\delta$-complete decision procedure. 
The main idea is to decide correctly whether slightly {\em relaxed} sentences over the reals 
are satisfiable or not. The following definitions are from \cite{DBLP:conf/lics/GaoAC12}.
\begin{definition}
A bounded quantifier is one of the following:
\[ \begin{array}{ll}
	\exists^{[a, b]} x  = \exists x : (a \le x \wedge x \le b) \\[1ex]
	\forall^{[a, b]} x  = \forall x : (a \le x \wedge x \le b)
\end{array} \]
\end{definition}

\begin{definition}
A bounded $\Sigma_1$ sentence is an expression of the form:
\begin{equation*} 
	\exists^{I_{1}} x_{1}, ..., \exists^{I_{1}} x_{n}: \psi(x_{1}, ..., x_{n})
\end{equation*}
where $I_{i} = [a_{i}, b_{i}]$ are intervals, $\psi(x_{1}, ..., x_{n})$ is a Boolean combination
of atomic formulas of the form $g(x_{1}, ..., x_{n})\, \texttt{op} \, 0$, where $g$ is a composition
of Type 2-computable functions and $\texttt{op} \in \{<, \le, >, \ge, =, \ne\}$.
\end{definition}
Any bounded $\Sigma_1$ sentence is equivalent to a $\Sigma_1$ sentence in which all the atoms are
of the form $f(x_{1}, ..., x_{n}) = 0$ (\ie, the only {\em \texttt{op}} needed is `=')
\cite{DBLP:conf/lics/GaoAC12}.
Essentially, Type 2-computable functions can be approximated arbitrarily well by finite
computations of a special kind of Turing machines (Type 2 machines); most of the `useful' functions 
over the reals (\eg, continuous functions) are Type 2-computable \cite{kobook}.

The notion of $\delta$-weakening \cite{DBLP:conf/lics/GaoAC12} of a bounded sentence is central
to $\delta$-satisfiability.
\begin{definition}
Let $\delta \in \mathbb{Q^+} \cup \{0\}$ be a constant and $\phi$ a bounded $\Sigma_1$-sentence 
in the standard form 
\begin{equation} \label{def:sigma1}
\phi = \exists^{I_{1}} x_{1}, ..., \exists^{I_{n}} x_{n}: \bigwedge^{m}_{i = 1}(\bigvee^{k_i}_{j = 1} f_{ij}(x_{1}, ..., x_{n}) = 0)
\end{equation}
where $f_{ij}(x_{1}, ..., x_{n}) = 0$ are atomic formulas. The $\delta$-{\em weakening} of $\phi$ 
is the formula:
\begin{equation*}
\phi^\delta = \exists^{I_{1}} x_{1}, ..., \exists^{I_{n}} x_{n}: \bigwedge^{m}_{i = 1}(\bigvee^{k_i}_{j = 1} |f_{ij}(x_1, ..., x_n)| \le \delta)
\end{equation*}
\end{definition}
Note that $\phi$ implies $\phi^\delta$, while the converse is obviously not true.
The bounded $\delta$-satisfiability problem asks for the following: given a sentence of 
the form (\ref{def:sigma1}) and $\delta \in \mathbb{Q^+}$, correctly decide between
\begin{itemize}
	\item \textbf{unsat}: $\phi$ is false,
	\item \textbf{$\delta$-true}: $\phi^\delta$ is true.
\end{itemize}
If the two cases overlap (\ie, $\phi$ is both false and $\delta$-satisfiable) then either decision 
can be returned, thereby causing a {\em false alarm}. Such a scenario reveals that the 
formula is {\em fragile} --- a small perturbation (\ie, a small $\delta$) can change the formula's
truth value. 
The dReal tool \cite{DBLP:conf/cade/GaoKC13} implements an algorithm for solving 
the $\delta$-satisfiability problem, \ie, a $\delta$-complete decision procedure. Basically, 
the algorithm combines a DPLL procedure (for handling the Boolean parts of the formula) with 
interval constraint propagation (for handling the real arithmetic atoms). 



bbal.tex
As an example of a hybrid system described above we consider the following scenario (Figure 1): 
the ball is launched from the initial point with initial speed $\upsilon_0$ and angle 
to the horizon $\alpha$. When the ball reaches the ground it reflects according to the laws 
of reflection.

\begin{figure}[ht!]
\centering
\includegraphics[width=90mm]{bouncing-ball.jpg}
\caption{Bouncing ball scenario}
\end{figure}

The ball trajectory is defined by the system of equations:
\begin{equation} \label{eq:Sx}
S_x(t) = S_{x_0} + \upsilon_0 t \cos{\alpha}
\end{equation}
\begin{equation} \label{eq:Sy}
S_y(t) = S_{y_0} + \upsilon_0 t \sin{\alpha} - \frac{g t^2}{2}
\end{equation}
where $S_x$ and $S_y$ are projections (on the axes $x$ and $y$, respectively) of the distance travelled by the ball starting at the point $(S_{x_0}; S_{y_0})$, $\upsilon_0$ is the initial speed, $\alpha$ is the angle to the horizon, $g$ is standard gravity and $t$ is time. In this example we assume that the initial point of the ball's centre is $(S_{x_0} = 0; S_{y_0} = 0)$

A hybrid system characterising the behaviour of the ball will consists of a single mode with the dynamics defined in (\ref{eq:Sx}) and (\ref{eq:Sy}) and the only jump to the same mode, during which the time $t$ is reset to $0$, the value of $S_x(t)$ before the jump is assigned to the initial value of $S_x(t)$ in the mode after the jump, and the speed of the ball is reduced by the coefficient $0.9$.

In order to comply with the description of the hybrid system with random continuous initial parameters, let $\alpha$ be uniformly distributed over the interval $[0, 0.5]$ with the probability density function $f_{A}(\alpha) = 2$. Let $\upsilon_0 = 20$ and $t \in [0, 3]$.

Let us consider a bounded probabilistic reachability problem: \textit{is the probability that 
the system reaches in 4 steps the box $D =  \{S_x \in [90, 91], S_y \in [1, 2]\}$ greater or equal 
to $0.062$}?

The Borel set of continuous parameters is formed as:
\begin{equation*}
\begin{split}
B =\{\alpha: \exists t_{0, q_0}, t_{1, q_0}, t_{2, q_0}, t_{3, q_0} \mathord{\in} [0, T]:\\
(S_{x0} = \upsilon_0 t_{0, q_0} \cos{\alpha}) \\
\wedge (S_{y0} = \upsilon_0 t_{0, q_0} \sin{\alpha} - \frac{gt_{0, q_0}^2}{2}) \wedge 
(S_{y0} = 0) \wedge (t_{0, q_0} > 0) \\
\wedge (S_{x1} = S_{x0} + 0.9 \upsilon_0 t_{1, q_0} \cos{\alpha}) \\
\wedge (S_{y1} = 0.9 \upsilon_0 t_{1, q_0} \sin{\alpha} - \frac{gt_{1, q_0}^2}{2}) \wedge 
(S_{y1} = 0) \wedge (t_{1, q_0} > 0) \\
\wedge (S_{x2} = S_{x1} + 0.9^2 \upsilon_0 t_{2, q_0} \cos{\alpha}) \\
\wedge (S_{y2} = 0.9^2 \upsilon_0 t_{2, q_0} \sin{\alpha} - \frac{gt_{2, q_0}^2}{2})
\wedge (S_{y2} = 0) \wedge (t_{2, q_0} > 0) \\
\wedge (S_{x3} = S_{x2} + 0.9^3 \upsilon_0 t_{3, q_0} \cos{\alpha}) \\
\wedge (S_{y3} = 0.9^3 \upsilon_0 t_{3, q_0} \sin{\alpha} - \frac{gt_{3, q_0}^2}{2}) \wedge \\
(S_{x3} \in [90, 91]) \wedge (S_{y3} \in [1, 2])\}
\end{split}
\end{equation*}

Functions $S_x(\alpha, t)$ and $S_y(\alpha, t)$ are invertible with respect to $\alpha$ on 
the domain $[0, 0.5]\times[0, 3]$. Then by applying Proposition \ref{prop:ivp} the formulated 
reachability problem can be encoded as a SMT formula and solved using dReal:
\begin{multline*} 
\exists \alpha \in B = [\alpha_{min}, \alpha_{max}]:\\
(F_A'(\alpha) = f_A(\alpha)) \wedge (F_A(\alpha_{min}) = 0) \wedge (F_A(\alpha) \ge 0.062) .
\end{multline*}
The actual dReal output on the formula above was as follows:
\begin{verbatim}
SAT with the following box:
 time_0:[0.03127620501613272,
                         0.03127624003371717]
 Fa_t:[0.0625524053632542,
                         0.06255247702735234]
 Fa_0:0
 a_0:[0.4292134593911809, 0.4292134922320248]
 t_10:[1.528734841443786, 1.528734859618606]
 t_20:[1.37586135729942, 1.375861377493665]
 t_30:[1.042516533544058, 1.042516555982107]
 t_0:[1.813824229846001, 1.813824261219404]
 t_1:[1.632441880862484, 1.63244191572182]
 t_2:[1.46919748799121, 1.469197526723806]
 t_3:[0.8312598060237518, 0.8312598490599685]
 t_00:[1.698594379660991, 1.698594396018329]
 a_t:[0.4604896994248981, 0.460489732265742]
\end{verbatim}

The probability of reaching the unsafe region by the system thus falls into the interval 
$[0.0625524053632542, 0.06255247702735234]$, that is, the interval associated to the final 
value of $F_A$ (corresponding to \verb#Fa_t# in the printout above). The obtained result was 
validated using a simple Monte Carlo method in MATLAB. The time variable was discretised over 
the interval $[0, 3]$ and the system was simulated using 10,000 uniformly random samples from 
the interval $[0, 0.5]$. The number of samples hitting the unsafe region was 620, thereby 
giving a {\em probabilistic} estimate of $\frac{620}{10,000}= 0.062$. In comparison, our method 
gives an interval of size $0.72\times 10^{-7}$ which is {\em guaranteed} to contain the actual 
probability.


deltareach.tex

The following definition of hybrid system is a slight variant of the standard one.
\begin{definition}\label{def:HS}
A hybrid system with {\em probabilistic} and {\em nondeterministic} initial parameters 
consists of the following components:
\begin{itemize}
\item $Q = \{q_0, ..., q_m\}$ a set of modes (discrete components of the system),
\item $D = D_0 \times \ldots \times D_p $ a domain of {\em discrete} random parameters, where 
	each $D_i$ is a finite set of reals,
\item $R = [r_1, s_1] \times \ldots \times [r_l, s_l] \subset \mathbb{R}^{l}$ a domain of 
	{\em continuous} random parameters,
\item $Z = [y_1, z_1] \times \ldots \times [y_o, z_o] \subset \mathbb{R}^{o}$ a domain of 
	{\em nondeterministic} parameters,
\item $X = [u_1, v_1] \times \ldots \times [u_n, v_n] \times [0, T] \subset \mathbb{R}^{n+1}$ a domain 
	of continuous variables,

\item $S = Q \times X$ is the {\em hybrid state space} of the system,
\item $\Lambda = D \times R \times Z$ is the {\em parameter space} of the system,
\item $U \subseteq S$ an {\em unsafe} region of the state space,
\end{itemize}
and predicates (or relations)
\begin{itemize}
	\item $\text{flow}_{q}(\lambda, \textbf{x}^{0}, \textbf{x}^{t})$ 
	mapping the parameter $\lambda \in \Lambda$ and the continuous state 
	$\textbf{x}^{0}$ at time 0 to state $\textbf{x}^{t}$ at time point $t \in [0, T]$ in mode $q$
	\item $\text{init}_q(\textbf{x}^{0})$ indicating that $s = (q, \textbf{x}^{0})$ belongs 
	to the set of initial states,
	\item $\text{jump}_{q \rightarrow q'}(\lambda, \textbf{x}^{t}, \textbf{x}^{0})$ indicating 
	that the system with parameter $\lambda\mathord{\in}\Lambda$ can make a transition from 
	mode $q$, upon reaching 
	the jump condition in continuous state $\textbf{x}^{t}$ at time $t \in [0, T]$, 
	to mode $q^\prime$ and setting the continuous state to $\textbf{x}^{0}$,
	\item $\text{unsafe}_{q}(\textbf{x}^{t})$ indicating that $s = (q, \textbf{x}^{t}) \in U$
\end{itemize}
For all $q\in Q$ the sets defined by $\text{flow}_q, \text{init}_q, \text{jump}_q$,
and $\text{unsafe}_q$ are Borel; $\text{flow}_q$ and $\text{jump}_q$ are restricted 
to be functions of $(\lambda,\textbf{x}^0)$ and $(\lambda,\textbf{x}^t)$, respectively.
\end{definition}
The parameters in $\Lambda$ are assigned in the initial mode and remain 
unchanged throughout the system's evolution. Also, the
Borel assumption for the sets defined by the predicates is a theoretical requirement for well-definedness
of probabilities, and in practice it is easily satisfied. 
%(We recall that Borel sets are obtained through countable union and complement of open sets.)
The continuous dynamics of the system is defined in each flow, and it can either be presented 
as a system of Lipschitz-continuous ODEs or explicitly. In this paper we focus on hybrid 
systems for which in each mode only one jump is allowed to take place (of course the model may have 
multiple jumps, but only one jump should be enabled at any time).
Given an initial value of the parameters, the semantics of a hybrid system can be informally 
thought as piece-wise continuous.
%: in each mode the system dynamics follows a continuous flow, 
%while a mode change (as result of a jump) might also involve a discontinuous change in the 
%variables $\textbf{x}$. 
%The evolution of a hybrid system is obtained by iteratively composing the relations {\em flow} and 
%{\em jump}, starting with {\em init} (more details can be found in \cite{DBLP:conf/hybrid/AlurCHH92}).
(More details about the formal semantics can be found in \cite{DBLP:conf/hybrid/AlurCHH92}.)

Bounded reachability asks whether the system reaches the unsafe region after $k \in \mathbb{N}$ 
discrete transitions.
\begin{definition}\cite{DBLP:journals/corr/GaoKCC14}\label{def:probreach}
The bounded $k$-step reachability property for hybrid systems with initial parameters
is the bounded $\Sigma_1$ sentence $ \exists \lambda\in\Lambda\  \psi (\lambda)$,
where
\begin{equation} \label{eq:reachability}
\begin{split}
\psi (\lambda) = \exists \textbf{x}^{0}_{0,q_0}, \exists \textbf{x}^t_{0,q_0}, ... , \exists \textbf{x}^{0}_{0,q_m}, \exists \textbf{x}^t_{0,q_m} , ... , \exists \textbf{x}^{0}_{k,q_m}, \exists \textbf{x}^t_{k,q_m}:\\
(\bigvee_{q \in Q} (\text{init}_{q}(\textbf{x}^{0}_{0,q}) \wedge \text{flow}_{q}(\lambda, \textbf{x}^{0}_{0,q},\textbf{x}^t_{0,q}))) \\
\wedge (\bigwedge_{i=0}^{k-1} (\bigvee_{q,q' \in Q}(\text{jump}_{q \rightarrow q'}(\lambda, \textbf{x}^t_{i,q}, \textbf{x}^{0}_{i+1,q'}) \\
\wedge (\text{flow}_{q'}(\lambda, \textbf{x}^{0}_{i+1,q'}, \textbf{x}^t_{i+1,q'}))) \wedge (\bigvee_{q \in Q} \text{unsafe}_{q}(\textbf{x}^t_{k,q}))))
\end{split}
\end{equation}
\end{definition}
%%%
\hide{
\begin{definition} \cite{DBLP:journals/corr/GaoKCC14}\label{def:reachability}
Bounded $k$-step $\delta$-reachability in hybrid systems can be encoded as a bounded $\Sigma_1$-sentence
\begin{equation} \label{eq:reachability}
\begin{split}
\exists \textbf{x}^{0}_{0,q_0}, \exists \textbf{x}^{t}_{0,q_0}, ... , \exists \textbf{x}^{0}_{0,q_m}, \exists \textbf{x}^{t}_{0,q_m} , ... , \exists \textbf{x}^{0}_{k,q_m}, \exists \textbf{x}^{t}_{k,q_m}:\\
(\bigvee_{q \in Q} (\text{init}_{q}(\textbf{x}^{0}_{0,q}) \wedge \text{flow}_{q}(\textbf{x}^{0}_{0,q},\textbf{x}^{t}_{0,q}))) \\
\wedge (\bigwedge_{i=0}^{k-1} (\bigvee_{q,q' \in Q}(\text{jump}_{q \rightarrow q'}(\textbf{x}^{t}_{i,q}, \textbf{x}^{0}_{i+1,q'}) \\
\wedge (\text{flow}_{q'}(\textbf{x}^{0}_{i+1,q'}, \textbf{x}^{t}_{i+1,q'}))) \wedge (\bigvee_{q \in Q} \text{unsafe}_{q}(\textbf{x}^{t}_{k,q}))))
\end{split}
\end{equation}
where $\textbf{x}^{0}_{i,q}$ and $\textbf{x}^{t}_{i,q}$ represent the continuous state at time point $t$ in the mode $q$ at the depth $i$, and $q'$ is a successor mode.
\end{definition}
}
\hide{
\begin{remark}
In this paper we are using a shorter version $\exists \textbf{x}^{t}_{i,q_j}$ of a longer notation $\exists i \in [0, k], \exists j \in [0, m], \exists q_{j} \in Q, \exists t \in [0, T], \exists \textbf{x}_{i,q_j}(t)$.
\end{remark}
}
%%%%
%Intuitively, the formula $\psi$ above can be understood as follows: the first conjunction is asking 
%for a set of continuous variables which satisfy the initial condition in one of the modes and the 
%flow in that mode; the second conjunction is looking for a set of vectors which satisfy any $k$ 
%discrete jumps and flows in each successor mode defined by the jumps; the third conjunction is 
%verifying whether the state of the system (the mode and the set of continuous variables in the 
%mode after $k$ jumps) belongs to the unsafe region. Note that the previous definition asks for
%reachability in {\em exactly} $k$ steps. 
Informally, the formula $\exists \lambda \in \Lambda\;\;  \psi(\lambda)$ encodes the sentence 
``there exists a parameter vector for which starting from {\em init} and following {\em flow} and {\em jump},
the system reaches the unsafe region in $k$ steps''.
We obtain reachability {\em within} $k$ steps by forming a disjunction of formula (\ref{eq:reachability})
for all values from 1 to $k$.
%The dReach tool can generate automatically such formulae from a hybrid system description.
The bounded reachability problem can be solved using a $\delta$-complete decision 
procedure \cite{DBLP:conf/lics/GaoAC12}, which will {\em correctly} return one of the following answers: 
\begin{itemize}
	\item \textbf{unsat}: meaning that formula (\ref{eq:reachability}) is unsatisfiable 
	(the system never reaches the bad region $U$);
	\item \textbf{$\delta$-sat}: meaning that formula (\ref{eq:reachability}) is 
	$\delta$-satisfiable. In this case a witness, \ie, an assignment for all 
	the variables, is also returned.
\end{itemize}
With a $\delta$-complete decision procedure, an \textbf{unsat} answer can always be trusted, 
while a \textbf{$\delta$-sat} answer might in fact be a false alarm caused by the overapproximation.
(In Appendix \ref{apndx:delta-sat} we provide a short overview of $\delta$-satisfiability.)

We now associate a probability measure to the random parameters, and we consider the following problem: 
what is the probability that a hybrid system with initial parameters reaches the unsafe region in $k$ steps?
Note that hybrid systems with both random and nondeterministic parameters will 
feature a range of reachability probabilities (although not necessarily a full interval). 
\begin{definition}\label{def:probboundreach}
The {\em probabilistic bounded $k$-step reachability problem} for hybrid system with initial parameters is
to compute an interval $[a,b]$ such that:
\begin{equation}\label{eq:probabilistic-reachability}
	\forall z_0\in Z \quad \int_{B|_{z_0}} \,dP \  \in \  [a,b]
\end{equation}
where 
\begin{equation} \label{def:Borel-set}
B = \{\lambda\in \Lambda : \psi(\lambda)\}
\end{equation}
and formula $\psi(\lambda)$ is per Definition \ref{def:probreach}; $P$ is the probability measure 
associated with the random parameters; and $B|_{z_0}$ is the restriction of $B$ to $z_0$.
\end{definition}
Informally, $B$ is the set of the parameter values for which the system reaches the unsafe region 
in $k$ steps. 
\begin{proposition}\label{prop:BBorel}
The set $B$ defined by (\ref{def:Borel-set}) is Borel.
\end{proposition}
(Proofs can be found in Appendix \ref{apndx:proofs}.)
The proposition entails that for any choice of the nondeterministic parameters, the probability that 
the system reaches the unsafe region is well-defined, and thereby Definition \ref{def:probboundreach} is
well-posed. When consider hybrid systems with continuous random parameters only, 
Definition \ref{def:probboundreach} can be strengthened.
\begin{definition}\label{def:probboundreach_rv}
Given any $\epsilon \in \mathbb{Q} \cap (0,1]$, the {\em probabilistic bounded $k$-step reachability problem}
for hybrid systems with {\em random continuous} initial parameters and single initial state is to compute 
an interval $[a,b]$ of length up to $\epsilon$ such that:
\begin{equation}\label{eq:probabilistic-reachability_rv}
	\int_{B} \,dP \  \in \  [a,b]
\end{equation}
where 
\begin{equation} 
B = \{\lambda\in \Lambda : \psi(\lambda)\}
\end{equation}
and formula $\psi(\lambda)$ is per Definition \ref{def:probreach}; $P$ is the probability measure 
associated with the random parameters.
\end{definition}
Note that if only discrete random parameters are present it might not be possible to obtain 
an arbitrarily small enclosure.
Also, in Definition \ref{def:HS} we require all continuous domains to be bounded: 
this is a necessary condition for $\delta$-decidability of bounded reachability \cite{DBLP:conf/lics/GaoAC12}. 
However, we later show that it is still possible to reason about random parameters with 
unbounded domains, \eg, normally distributed. The key is that any probability density function
can be approximated arbitrarily well by a truncation on a large (but finite) interval. 


%If the random parameters are independent then 
%$P(\textbf{r}) = \Pi_{i=1}^{l} P_i(r_{i})$, where $P_i$ is the probability measure of 
%parameter $r_i$. The infinitesimal probability measure $dP(\textbf{r})$ is usually 
%given as a (probability) density function.

integrate.tex
In order to compute probabilistic reachability in hybrid systems with unbounded random parameter it is sufficient to solve bounded $\Sigma_1$ sentences of the following type:
\begin{equation} \label{eq:integration-problem}
\exists x \in [a, b] : \int_{a}^{x} f(t)\, dt \ge C
\end{equation}
for some constant $C\in [0,1]$.

\begin{proposition} \label{prop:ivp}
Let $f\mathord{:}[a,b] \rightarrow \mathbb{R}$ be a Lipschitz-continuous function. Then there 
exists one and only one continuous and differentiable function $F$ on $[a, b]$ 
such that $F(a) = 0$ and $\forall x \in [a, b]: \int_a^x f(t)\, dt = F(x)$
\end{proposition}

\begin{proof}
Let us formulate an initial value problem:
\begin{equation} \label{eq:IVP}
(F'(x) = f(x)) \wedge (F(a) = 0)
\end{equation}
where $f(x)$ is Lipschitz-continuous in $F(x)$ and continuous on $[a,b]$. 
By the Picard-Lindel\"{o}f theorem and the assumption $F(a)=0$, the formulated IVP 
has a unique solution on $[a, b]$ derived as following:
\begin{equation*} 
\forall x \in [a,b]: \int_a^x f(t)\, dt = F(x) - F(a) = F(x)
\end{equation*}
\end{proof}

Proposition \ref{prop:ivp} allows assessing the value of the integral by the value of the 
function which is obtained from the initial value problem (\ref{eq:IVP}). Therefore, 
the integration problem (\ref{eq:integration-problem}) is equivalent to the sentence
\begin{equation*} %\label{eq:encoding}
\exists x \in [a, b] : (F'(x) = f(x)) \wedge (F(x) \ge C) \wedge (F(a) = 0)
\end{equation*}
which can be directly solved by dReal as a $\delta$-satisfiability problem.


intro.tex

Cyber-physical systems integrates digital computing (the {\em cyber} part) with a {\em physical} 
environment or device, in order to enhance or enable new capabilities 
of physical systems. Hybrid systems are mathematical models that combine continuous dynamics 
and discrete control, and enjoy widespread use for modelling cyber-physical systems.
For example, Stateflow/Simulink\footnote{\url{www.mathworks.com/simulink}} is the {\em de facto} 
standard tool for model-based design of embedded systems, and its semantics can be given 
in terms of hybrid systems (\eg, \cite{Tiwari02}). Cyber-physical systems are
used in many safety-critical applications, where a malfunctioning can result in threats to, 
or even loss of, human life. For example, modern aircraft are flown more efficiently 
by a computer, while anti-lock brakes and stability control
contribute to safer cars. Again, electronic biomedical devices (\eg, digital infusion pumps) 
%that deliver drugs or nutrients to patients) 
offer superior flexibility and accuracy than traditional devices. Thus, verifying safety of 
cyber-physical systems, and thereby of hybrid systems, is an extremely important problem.

The state space of a hybrid system consists of a discrete component and of a continuous component. 
The fundamental {\em reachability} problem is to decide whether a hybrid system reaches an 
{\em unsafe} region of its state space (a subset of states indicating incorrect behaviour of 
the system). Unfortunately, this problem is undecidable even for hybrid systems with constant 
differential dynamics \cite{DBLP:conf/hybrid/AlurCHH92}. 
For timed automata, \ie, same constant differential dynamics across all the variables, 
the reachability problem is PSPACE-complete \cite{DBLP:conf/icalp/AlurD90}. Also, it has been 
recently shown that bounded-time reachability of rectangular automata with non-negative rates 
is decidable \cite{BrihayeDGORW11}.
However, hybrid systems arising from practical applications feature much richer dynamics, including
non-linear functions over the reals, \eg, trigonometric functions,
for which even simple questions are in general undecidable \cite{Richardson68}.
Furthermore, for many practical applications it is necessary to augment hybrid systems with 
stochastic behaviour. Stochastic systems arise naturally when modelling phenomena which are 
intrinsically probabilistic, \eg, soft errors in computing hardware. Also, stochastic systems 
can arise due to uncertainty in (deterministic) system components, its behaviours, and its 
environment. The reachability problem for stochastic hybrid systems asks what is the {\em probability} 
that the system reaches the unsafe region. (Note that for hybrid systems with both stochastic and 
non-deterministic behaviour the answer may be a range of probabilities.)
In this work we focus on {\em bounded} reachability, \ie, within a finite number of discrete 
transitions.

Since even standard reachability is undecidable, the problem must be modified if we want 
to solve it algorithmically. A possible solution is to relax it in a sound manner through
the notions of $\delta$-satisfiability and $\delta$-complete decision 
procedures \cite{DBLP:conf/lics/GaoAC12}. Such procedures sidestep undecidability by allowing a
`tuneable' precision in the answer provided. This is a necessary condition for decidability, 
and it motivates the notion 
of $\delta$-satisfiability for logical formulae over the reals \cite{DBLP:conf/lics/GaoAC12}. 
Using $\delta$-satisfiability, in this paper we introduce and study the notion of probabilistic 
$\delta$-reachability.
%Typically, $\delta$-complete
%procedures employ SMT solvers and rigorous numerical methods 
%to verify an overapproximation of the system behaviour, so that bugs will always be detected. 
%Although the over approximation can be tight (tuneable by an arbitrarily small rational parameter, 
%$\delta$), the decision procedure may produce false alarms. 
%It may report a bug when in reality 
%the system is safe. However, such a scenario indicates that the system is {\em fragile},
%thereby providing valuable information to the system 
%designer. 
%The dReal tool \cite{DBLP:conf/cade/GaoKC13} implements $\delta$-complete procedures for 
%first-order formulae that may include solutions of nonlinear ordinary differential equations (ODEs).

To summarise, in this paper:
\begin{itemize}
	\vspace{-1ex}
	\item we formulate the bounded $\delta$-reachability problem for hybrid systems 
		with continuous/discrete probabilistic and nondeterministic initial parameters;
	\item we develop an algorithm that combines validated integration and $\delta$-complete 
		procedures into for computing a {\em numerically guaranteed} enclosure for the 
		reachability probabilities. For models with continuous random 
		(but no nondeterministic) parameters, such enclosure can be made {\em arbitrarily} small;
		%mention ODEs
	\item we validate our algorithm against standard Monte Carlo probability estimation 
		on a number of case studies.
\end{itemize}

\paragraph{{\bf Related Work.}}
The SiSAT tool \cite{Sisat} solves probabilistic bounded 
reachability by returning answers guaranteed to be numerically accurate.
However, SiSAT does not currently support continuous random parameters, while instead our
tool does so (also with unbounded domains, \eg, normal random variables). A very recent extension 
of SiSAT supports continuous nondeterminism, but the technique is based on statistical model checking 
and therefore can only provide statistical guarantees \cite{SMCnd}, while we give numerical and formal
guarantees.
In \cite{Enszer10verifiedsolution} the authors present a technique for computing p-boxes using validated ODE
integration. However, the technique is restricted to ODE systems and finite-support random parameters, 
while we handle hybrid systems and infinite-support random parameters. Moreover, it is 
not clear what guarantees are given for models containing only continuous and/or discrete random parameters: 
the size of the computed p-box might be quite large. In contrast, for continuous random parameters we can 
compute an arbitrarily small interval containing the exact reachability probability (see Proposition \ref{prop:cn_nd}).

UPPAAL \cite{UPPAAL} is an 
extremely powerful model checker for timed automata, and it has been recently extended to support 
(dynamic) networks of stochastic timed automata via UPPAAL SMC \cite{UPPAALSMCtut}. However,
UPPAAL SMC utilises a statistical model checking approach for reasoning
about probabilities. PRISM \cite{KNP11} is a state-of-the-art model checker for a variety
of discrete-state stochastic systems, but with respect to real-time systems it is limited to
probabilistic timed automata. The tool {FAUST}$^2$ \cite{faust2} utilises abstraction
techniques to verify nondeterministic continuous-state Markov models, although currently
for discrete-time models only. ProHVer computes an upper bound for the maximal reachability 
probability \cite{ZhangSRHH10}, and handles continuous random parameters via discrete overapproximation 
only \cite{FraenzleHHWZ11}. We instead provide an enclosure ({\em both} upper and lower bounds) of the 
whole range of probabilities (for models with nondeterministic continuous parameters); in the case of 
continuous random parameters our enclosure can be arbitrarily tight (see Proposition \ref{prop:cn_nd}).
In \cite{Abate2010624} the authors introduce a technique for computing bounds on reachability probababilities 
for stochastic hybrid systems, using abstraction by discrete-time Markov chains.
The technique is further extended to full LTL and nondeterminism \cite{TkachevA13}.
In \cite{RamponiCSL10} the authors give model checking algorithms for PCTL formulae over continuous-time 
stochastic hybrid systems. However, in \cite{Abate2010624,TkachevA13,RamponiCSL10} continuous state space 
is handled through finite discretisation and approximated numerical solutions are provided for the experiments. 
We instead consider continuous time and space, and give full mathematical/numeric guarantees.

With respect to $\delta$-satisfiability, in \cite{Novak1992} the authors
introduced and studied the complexity of a {\em relaxed} version of the verification problem, \ie, verifying whether
a given candidate is close to a problem solution. The (strong) verification problem is 
undecidable in general, so the authors relax it by introducing a ``safety zone'' in
which either answer is deemed correct --- this is $\delta$-satisfiability. 
%The authors derived complexity bounds for solving the relaxed verification problem for operators 
%defined over vector spaces and computed via linear functionals.
%Novak and Wo\'{z}niakowski worked in the context of information-based complexity \cite{CandI} 
%and derived general complexity bounds for solving the relaxed verification problem for operators 
%defined over vector spaces and computed via linear functionals (thereby including differentiation 
%and integration).

Finally, in our work we use verified integration techniques (for an overview see, \eg, 
\cite{Petras:VerifiedNI} and references therein). Integration methods in the literature work 
with integrands in explicit form, \ie, one must provide the actual mathematical expression 
for the integrand. Our approach is more general because:
a) the integrand is given as a function of the numerical solution of possibly nonlinear ODEs; 
b) it considers hybrid dynamics.  Our algorithm carries over the guarantees provided by 
$\delta$-complete procedures for aspects a) and b) to the verified computation of a 
multi-dimensional integral over a possibly unbounded domain. 


main.tex
\documentclass[envcountsame]{llncs}

\usepackage{amssymb}
\usepackage{url}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{fancyvrb}

\usepackage[boxed,ruled,vlined]{algorithm2e}
%\usepackage{relsize}
%\usepackage{makeidx} 

\usepackage{multirow}

% allows for temporary adjustment of side margins
\usepackage{chngpage}

%\pagestyle{headings} 

\newcommand{\hide}[1]{}
\newcommand{\eg}{{\em e.g.}}
\newcommand{\ie}{{\em i.e.}}
\newcommand{\preach}{{\sffamily ProbReach}}

\pagestyle{plain}


\begin{document}

\title{Probabilistic bounded reachability for hybrid systems with continuous nondeterministic 
and probabilistic parameters}

\author{Fedor Shmarov \and Paolo Zuliani}
\institute{School of Computing Science, Newcastle University, Newcastle upon Tyne, UK \\
\email{\{f.shmarov, paolo.zuliani\}@ncl.ac.uk}}

\maketitle

\begin{abstract}
%Verification of cyber-physical systems is a difficult, yet extremely important, problem.
%Hybrid systems offer a modelling framework in which to perform formal verification of
%cyber-physical systems. 
We develop an algorithm for computing bounded reachability
probability for hybrid systems, i.e., the probability that the system reaches an unsafe
region within a finite number of discrete transitions. In particular, we focus on hybrid
systems with continuous dynamics given by solutions of nonlinear ordinary 
differential equations (with possibly nondeterministic initial conditions and parameters), and 
probabilistic behaviour given by initial parameters distributed as continuous (with 
possibly infinite support) and discrete random variables. 
Our approach is to define an appropriate relaxation of the (undecidable) reachability problem, 
so that it can be solved by $\delta$-complete decision procedures. In particular, for systems 
with continuous random parameters only, we develop a validated integration procedure which 
computes an arbitrarily small interval that is guaranteed to contain the reachability probability. 
In the more general case of systems with both nondeterministic and probabilistic parameters, our 
procedure computes a guaranteed enclosure for the range of reachability probabilities. 
We have applied our approach to a number of nonlinear hybrid models and validated the results 
by comparison with Monte Carlo simulation.

%Moreover, we discuss certain types of random variables with unbounded support and show 
%that the bounded $\delta$-reachability problem can still be solved by using an appropriate 
%$\delta$-complete decision procedure. 

\end{abstract}

\section{Introduction}
\input{intro}

\section{Probabilistic Bounded $\delta$-Reachability}
\input{deltareach}

%\input{probdeltareach}

%\section{Systems with one random parameter}
%\input{onerand}

\input{probreachimpl}

%\section{Integration problem via differentiation}
%\input{integrate}

%%\input{multivar}

\input{algorithm}

\input{results}

\vspace{-1ex}
\section{Conclusions and Future Work}
We have given a formal definition of the bounded probabilistic $\delta$-reachability problem
for hybrid systems with continuous random and nondeterministic initial parameters. We have combined
validated integration with $\delta$-complete decision procedures for solving the 
probabilistic $\delta$-reachability problem. Our technique computes a {\em numerically guaranteed}
enclosure for the probabilities that the system reaches the unsafe region in a finite number of 
discrete transitions. For systems with continuous random (but no nondeterministic) parameters,
such enclosure can be made arbitrarily small.
We have implemented our technique in the open source tool \preach\  and have applied
it to a number of case studies featuring highly nonlinear ODEs, unbounded continuous random 
parameters and nondeterministic parameters. We have validated our results against Monte Carlo
simulation, and the comparison supports the correctness of our approach. 
\hide{
Also, our technique can produce faster and accurate results, returning intervals 
of size $10^{-9}$ even when computing very small probabilities (see, \eg,
the results for the bouncing ball model). This is a significant advantage with respect 
to Monte Carlo verification techniques such as statistical model checking
\cite{YounesS06}. It is well-known that Monte Carlo methods suffer from the rare-event 
problem: to estimate reliably very small probabilities, 
extremely large sample sizes (\ie, system simulations) are needed \cite{RRbook,Z:hscc12}.
}

Our work shows that it is possible to verify bounded reachability for hybrid systems 
featuring continuous random and nondeterministic parameters with the same level of accuracy as
for finite-state stochastic systems. Of course, more work needs to be done in terms of improving both 
the tool engineering and the theory. With respect to the former, a more efficient parallel strategy
needs to be implemented, and more experiments need to be performed to assess better the tool scalability. 
For the theory, in the future we plan to tackle a larger class of hybrid systems, which in particular
include state-dependent probabilistic jumps and continuous probabilistic dynamics (stochastic 
differential equations).
\hide{
A first extension is to allow probabilistic jumps in the system (discrete) dynamics. Furthermore, we plan 
to tackle continuous random dynamics such as, \eg, stochastic differential equations (SDEs). These are 
used, \eg, for modelling white noise in control systems, and are therefore important for modelling
cyber-physical systems. Thus, solving bounded $\delta$-reachability in general stochastic 
hybrid systems requires building a $\delta$-complete decision procedure for SDEs, which would in
turn require a validated SDE solver --- this is a very challenging problem.
Another interesting direction of research would be combining semi-exhaustive methods to search
the state space (\eg, quasi-Monte Carlo methods) with our approach. As the number of boxes to analyse 
grows exponentially with the number of initial parameters, any technique that would accelerate 
convergence is naturally of great interest.
}


\hide{
We believe that our technique can be fruitfully utilised for model 
checking \cite{DBLP:conf/lop/ClarkeE81} actual implementations of probabilistic algorithms. 
In particular, by adopting a similar approach as with bounded model checking for 
C code \cite{ckl2004}, temporal logic properties over programs would be transformed into
SMT formulae by `unrolling' the program source. Of course, these programs should satisfy suitable 
conditions such as, \eg, having finite loops only. We think that our approach would be useful for
verifying embedded software implementing nonlinear control laws in (continuous) environments
described by complex ODEs.
}
\hide{\vspace{-2ex}
\paragraph{Acknowledgement}
This work has been supported by award N00014-13-1-0090 of the US Office of Naval Research.}


\bibliographystyle{splncs03}
\vspace{-3ex}
\bibliography{../refs}

\newpage
\appendix
\section*{Appendix}
\input{appendix}

\end{document}



multivar.tex
\section{Multiple random variables}
In this section we discuss the enhancements to the proposed algorithm for solving probabilistic bounded reachability in hybrid systems with multiple independent random continuous initial parameters. Essentially, the structure of the algorithm will remain the same but the error estimation and box division strategy will be changed.

\subsection{Error estimation}
In order to calculate the probability of reaching an unsafe region by a hybrid system with multiple independent continuous random initial parameters, it is sufficient to obtain a Cartesian product of partitions over each random variable and verify each {\em box} from the Cartesian product.
\begin{proposition}
Given a hybrid system with $l$ independent continuous random initial parameters and a desired size of the probability interval $\epsilon_{prod} \in \mathbb{Q}^{+}$ each random variable should be integrated with precision $\epsilon$ satisfying the formula below:
\begin{equation} \label{eq:multivarmaininequalityspecialcase}
\begin{split}
\epsilon_{prod} \le {{l}\choose{1}}\epsilon + {{l}\choose{2}}\epsilon^{2} + ... + {{l}\choose{i}}\epsilon^{i} + ... + {{l}\choose{l-1}}\epsilon^{l-1}  + {{l}\choose{l}}\epsilon^{l}
\end{split}
\end{equation}
where ${{l}\choose{i}}$ - binomial coefficient.
\end{proposition}
\begin{proof}
If a hybrid system has $l$ independent initial random parameters then the probabilistic reachability can be computed as:
\begin{equation} \label{eq:mulvar}
\int_{B} \prod_{i = 1}^{l} dP(r_{i})
\end{equation}
where $dP_{i}$ is a probability measure of the $i$-th random parameter and $B$ is a Borel set such that $B = [a_{1}, b_{1}] \times ... \times [a_{l}, b_{l}]$ containing all the values of random parameters for which the hybrid system reaches the unsafe region.

We can employ the technique described above for computing the formula (\ref{eq:mulvar}). The formula will be presented in the form:
\begin{equation} \label{eq:mulvarind}
\int_{\Omega} I_{B}(r_{1}, ..., r_{l}) \prod_{i = 1}^{l} dP(r_{i})
\end{equation}
where $I_{B}(r_{1}, ..., r_{l})$ is the indicator function which is equal to 1 if the system reaches the unsafe region and 0 otherwise.

Using the definition of the indicator function the formula (\ref{eq:mulvarind}) is equivalent to
\begin{equation}
\int_{\Omega} \prod_{i = 1}^{l} dP(r_{i})
\end{equation}
for all $r_{i}$ for which the system reaches the unsafe region. The value above can be calculated as a product of the integrals.

\begin{equation}
\int_{\Omega} \prod_{i = 1}^{l} dP(r_{i}) = \prod_{i = 1}^{l} \int_{a_{i}}^{b_{i}} dP(r_{i}) = \prod_{i = 1}^{l} I_{i}
\end{equation}

Now it is desired for the product of the integral to be contained in the interval of defined length $\epsilon_{prod}$. Taking into account that we can obtain an interval of length $\epsilon_{i}$ containing its {\em exact} value of the $i$-th integral its sufficient to demonstrate how the values $\epsilon_{i}$ should be chosen so that the product of integrals is contained within the interval of length $\epsilon_{prod}$.

According to the rules of interval arithmetics product of the integrals is contained in the interval:
\begin{equation}
[\prod_{i = 1}^{l} I_{i}^{est}, \prod_{i = 1}^{l} I_{i}^{est} + \epsilon_{prod}] = [I_{1}^{est}, I_{1}^{est} + \epsilon_{1}] \cdot ... \cdot [I_{l}^{est}, I_{l}^{est} + \epsilon_{l}] = [\prod_{i = 1}^{l} I_{i}^{est}, \prod_{i = 1}^{l} (I_{i}^{est} + \epsilon_{i})]
\end{equation}

Therefore, given $\epsilon_{prod}$ the values of $\epsilon_{i}$ should be chosen in such a way that the following inequality holds:
\begin{equation} \label{eq:mulvarinequality}
\prod_{i = 1}^{l} I_{i}^{est} + \epsilon_{prod} \le \prod_{i = 1}^{l} (I_{i}^{est} + \epsilon_{i})
\end{equation}

\iffalse
Let $R$ denote a set $\{I_{1}^{est}, ..., I_{n}^{est}\}$ and $S = \{\epsilon_{1}, ..., \epsilon_{n}\}$. Then the right-hand side of the formula (\ref{eq:mulvarinequality}) can be modified:
\begin{equation}
\begin{split}
\sum ({{R}\choose{0}}{{S}\choose{l}}) + \sum ({{R}\choose{1}}{{S}\choose{l - 1}}) + ... + \sum ({{R}\choose{i}}{{S}\choose{l - i}}) + \\
... + \sum ({{R}\choose{l - 1}}{{S}\choose{1}}) + \sum ({{R}\choose{l}}{{S}\choose{0}})
\end{split}
\end{equation}
where $\sum ({{R}\choose{i}}{{S}\choose{l - i}})$ is a sum of products of all the combinations of length $i$ from the set $R$  multiplied by the combinations of length $l - i$ from set $S$ chosen in such a way that all the indexes in each of the obtained products do not repeat. 

Taking into account that $I_{i}^{est} \in [0, 1]$ the following holds:
\begin{equation} \label{eq:multivarmaininequality}
\begin{split}
\epsilon_{prod} \le \sum {{S}\choose{l}} + \sum {{S}\choose{l - 1}} + ... + \sum {{S}\choose{l - i}} + \\
... + \sum {{S}\choose{1}} + \sum {{S}\choose{0}}
\end{split}
\end{equation}
\fi

Therefore, choosing $\epsilon_{i}$ in such a way that (\ref{eq:mulvarinequality}) holds will guarantee that the {\em exact} value of the product of $l$ integrals is contained in the interval of size $\epsilon_{prod}$. In a special case when all $\epsilon_{i}$ are equal to $\epsilon$ the formula (\ref{eq:mulvarinequality}) can be simplified:
\begin{equation} \label{eq:multivarmaininequalityspecialcase}
\begin{split}
\epsilon_{prod} \le {{l}\choose{1}}\epsilon + {{l}\choose{2}}\epsilon^{2} + ... + {{l}\choose{i}}\epsilon^{i} + ... + {{l}\choose{l-1}}\epsilon^{l-1}  + {{l}\choose{l}}\epsilon^{l}
\end{split}
\end{equation}
where ${{l}\choose{i}}$ is a binomial coefficient
\end{proof}

\subsection{Division strategy}
Given a box $D = [a_{1}, b_{1}] \times \cdot \times [a_{n}, b_{n}]$ we can partition it into $2_{n}$ boxes of the equal sizes in such a way that each interval $[a_{i}, b_{i}]$ is altered. This procedure should be applied to every box for which decision procedure defined above returns \textbf{$\delta$-sat} for both formulas until the desired precision is reached.

\begin{remark}
Essentially, any division strategy can be applied as long as the size of each interval forming the box is reduced.
\end{remark}


onerand.tex
In the scope of this paper we consider hybrid systems with one random initial parameter, and 
we propose an algorithm solving the probabilistic $\delta$-reachability problem for such 
hybrid systems. In a hybrid system with random parameters, the random variable $r$ is defined 
over a bounded interval. However, we will demonstrate that we can correctly evaluate 
probabilistic reachability problems also for certain types of unbounded random variables.

\begin{proposition} \label{prop:inf}
The probabilistic reachability problem for hybrid systems with a single random parameter 
can be correctly solved by a $\delta$-complete decision procedure if:
\begin{itemize}
	\item the domain of the random parameter is bounded from one side,
	\item the domain of the random parameter is unbounded and its probability density 
		function is symmetric.
\end{itemize}
\end{proposition}
\begin{proof}
Let's consider the first case when the domain $\Omega$ of the continuous random variable is bounded from one side. The problem when set $B \subseteq \Omega$ is bounded is trivial and probabilistic reachability problem can be solved correctly by the $\delta$-complete procedure.

If $B$ is unbounded then its complement $\Omega/B$ is bounded (because $\Omega$ is bounded from one side and $B \subseteq \Omega$). Set $\Omega/B$ can be constructed by inverting the relation $\text{unsafe}_{q}(\textbf{r},\textbf{x})$. Taking into account that
\begin{equation} \label{eq:total-probability}
\int_{B} \,dP(r) + \int_{\Omega/B} \,dP(r) = \int_{\Omega} \,dP(r) = 1
\end{equation}
probabilistic reachability problem (\ref{eq:probabilistic-reachability}) is equivalent to
\begin{equation} \label{eq:probabilistic-reachability-one-side}
\exists B_1 \subseteq  \Omega/B: 1 - \int_{B_1} \,dP(r) < C
\end{equation}

Therefore, solving two problems (\ref{eq:probabilistic-reachability}) and (\ref{eq:probabilistic-reachability-one-side}) at the same time by the $\delta$-complete decision procedure guaranties that one of them will terminate with the correct result. Knowing one, we can derive another from the equation (\ref{eq:total-probability}).

Let's consider the second case when $\Omega$ is unbounded and the probability density function is symmetric. We can find two disjoint bounded from one side sets $\Omega^{-}, \Omega^{+} \subset \Omega$ such that $\Omega^{-} \cup \Omega^{+} = \Omega$ and $\int_{\Omega^{-}} \,dP(r) = \int_{\Omega^{+}} \,dP(r) = \frac{\int_{\Omega} \,dP(r)}{2} = 0.5$.

Let $B \subseteq \Omega$. We can find two disjoint sets $B^{-}, B^{+} \subset B$ such that $B^{-} \cup B^{+} = B$, $B^{-} \subseteq \Omega^{-}$ and $B^{+} \subseteq \Omega^{+}$. Then problem (\ref{eq:probabilistic-reachability}) is equivalent to:
\begin{equation} \label{eq:probabilistic-reachability-symmetric}
\exists B_{1}^{-} \subseteq  B^{-}, \exists B_{1}^{+} \subseteq  B^{+} : \int_{B_{1}^{-}} \,dP(r) + \int_{B_{1}^{+}} \,dP(r) \ge C
\end{equation}
where for each of the subsets $B^{-}$ and $B^{+}$:
\begin{equation}
\begin{split}
\int_{B^{+}} \,dP(r) + \int_{\Omega^{+}/B^{+}} \,dP(r)  = 0.5 \\
\int_{B^{-}} \,dP(r) + \int_{\Omega^{-}/B^{-}} \,dP(r)  = 0.5
\end{split}
\end{equation}
If either of subsets $B^{-}$ and $B^{+}$ is unbounded then the technique described above should 
be applied.
\end{proof}

In general, the set $B$ obtained in (\ref{def:Borel-set}) will be of the form:
\begin{equation*} %\label{set-B-one-parameter}
B = \bigcup_{i=1}^{w} B_{i}
\end{equation*}
for some $w \in \mathbb{N}$,
where $B_{i} = [\underline{r}_{i}, \overline{r}_{i}]$ and 
$\forall i \ne j \in \{1,\ldots, w\}:B_{i} \cap B_{j} = \emptyset$. However, if the functions 
representing the dynamics of the system are invertible with respect to the random parameter, 
then set $B$ is formed by a single interval. This allows us to write the computation of integral 
(\ref{eq:probabilistic-reachability}) as a bounded $\Sigma_1$ sentence 
\begin{equation*} %\label{eq:probability-one-variable}
\exists r \in [a, b]: \int_{a}^{r} f_R(x)\, dx \ge C
\end{equation*}
where $f_R$ is the probability density function of the random parameter and $C \in [0,1]$ is a
constant.

The probabilistic bounded reachability problem in hybrid systems with one random initial parameter 
and with invertible dynamics with respect to the random variable can be solved using a validated 
ODE solver that allows encoding integrals as IVPs (initial value problem). We aim for this 
formalisation as that is the way in which dReal handle ODE dynamics.




probdeltareach.tex
%We now introduce $\delta$-reachability in a probabilistic context.

\begin{definition}\label{def:pHS}
A hybrid system with {\em continuous random initial parameters} consists of $Q$, $X$ and $U$ as in
Definition \ref{def:HS}, and the following components:
\begin{itemize}
%\item $Q = \{q_0, ..., q_m\}$ a set of modes (discrete components of the system),
%\item $X = [u_1, v_1] \times ... \times [u_n, v_n] \times [0, T] \subset \mathbb{R}^{n+1}$ a domain 
	%of continuous variables,
\item $R = [a_1, b_1] \times ... \times [a_l, b_l] \subset \mathbb{R}^{l}$ a domain of {\em random}
	parameters,
\item $S = Q \times R \times X$ is the state space of the system,
%\item $U \subseteq S$ an {\em unsafe} region of the state space,
\end{itemize}
predicates $\text{flow}_{q}(\textbf{r}, \textbf{x}^{0}, \textbf{x}^{t})$,
$\text{init}_q(\textbf{r}, \textbf{x}^{0})$, 
$\text{jump}_{q \rightarrow q'}(\textbf{r}, \textbf{x}^{t}, \textbf{x}^{0})$ and
$\text{unsafe}_{q}(\textbf{r}, \textbf{x}^{t})$.
For all $q\in Q$ the sets defined by $\text{flow}_q, \text{init}_q, \text{jump}_q$,
and $\text{unsafe}_q$ are Borel.
\end{definition}
The random parameters $\mathbf{r}$ are assigned in the initial mode and remain 
unchanged throughout the system's evolution. Also, the
Borel assumption for the sets defined by the predicates is a theoretical requirement for well-definedness
of probabilities, and in practice it is easily satisfied. (We recall that Borel sets are obtained
through countable union and complement of open sets.)

We now modify the bounded $k$ step $\delta$-reachability problem as follows.
\begin{definition}\label{def:probreach}
The bounded $k$ step $\delta$-reachability property for hybrid systems with random initial parameters 
is the bounded $\Sigma_1$ sentence:
\begin{equation*} 
\exists \textbf{r}\   \psi (\textbf{r})
\end{equation*}
where
\begin{equation*} 
\begin{split}
\psi (\textbf{r}) = \exists \textbf{x}^{0}_{0,q_0}, \exists \textbf{x}^t_{0,q_0}, ... , \exists \textbf{x}^{0}_{0,q_m}, \exists \textbf{x}^t_{0,q_m} , ... , \exists \textbf{x}^{0}_{k,q_m}, \exists \textbf{x}^t_{k,q_m}:\\
(\bigvee_{q \in Q} (\text{init}_{q}(\textbf{r}, \textbf{x}^{0}_{0,q}) \wedge \text{flow}_{q}(\textbf{r}, \textbf{x}^{0}_{0,q},\textbf{x}^t_{0,q}))) \\
\wedge (\bigwedge_{i=0}^{k-1} (\bigvee_{q,q' \in Q}(\text{jump}_{q \rightarrow q'}(\textbf{r}, \textbf{x}^t_{i,q}, \textbf{x}^{0}_{i+1,q'}) \\
\wedge (\text{flow}_{q'}(\textbf{r}, \textbf{x}^{0}_{i+1,q'}, \textbf{x}^t_{i+1,q'}))) \wedge (\bigvee_{q \in Q} \text{unsafe}_{q}(\textbf{r}, \textbf{x}^t_{k,q}))))
\end{split}
\end{equation*}

\end{definition}

If we associate a probability measure to the random parameters, then we can assess {\em quantitative} 
system properties (such as the probability that the system reaches the unsafe region).
In particular, we consider the following problem: 
\textit{what is the probability that a hybrid system with random initial parameters reaches the 
unsafe region in $k$ steps?}
To calculate such probability, we first need to find the set 
that contains all the values of the random variables satisfying the reachability property. 
Such a set is constructed as follows:
\begin{equation} \label{def:Borel-set}
B = \{\textbf{r} : \psi(\textbf{r})\}
\end{equation}
where formula $\psi(\textbf{r})$ is per Definition \ref{def:probreach}.
\begin{proposition}\label{prop:BBorel}
The set $B$ defined by (\ref{def:Borel-set}) is Borel.
\end{proposition}
(Proofs can be found in the Appendix.)

The proposition entails that the probability that the system reaches the unsafe region is well-defined. 
Such a probability is computed by integrating the probability measure of  the random variables over 
set $B$. In particular, we need to compute the following integral
\begin{equation} \label{eq:probabilistic-reachability}
\int_{B} \,dP(\textbf{r})\quad  (= \int_{\Omega} I_{U}(\textbf{r})dP(\textbf{r}))
\end{equation}
where $B$ is the set (\ref{def:Borel-set}), $P(\textbf{r})$ is a probability measure over 
the random parameters $\textbf{r}$, $\Omega$ is the support of the random parameters, and 
$I_{U}$ is the indicator function
\begin{equation*}
        I_{U}(\textbf{r}) = 
        \begin{cases}
                1 \quad \text{if the system with parameter $\textbf{r}$ reaches $U$ in $k$-steps} \\
                0  \quad\text{otherwise.}
        \end{cases}
\end{equation*}
If the random parameters are independent then 
$P(\textbf{r}) = \Pi_{i=1}^{l} P_i(r_{i})$, where $P_i$ is the probability measure of 
parameter $r_i$. The infinitesimal probability measure $dP(\textbf{r})$ is usually 
given as a (probability) density function.


probreachimpl.tex
%\section{Probabilistic $\delta$-Reachability}
\section{Validated Integration Procedure}\label{sec:vip}
We now present the first part of our $\delta$-complete procedure for calculating the $k$-step
reachability probability (\ref{eq:probabilistic-reachability_rv}).
The algorithm consists of a validated integration procedure and a decision procedure used 
for computing the set $B$ of Definition \ref{def:probboundreach_rv}. For clarity, we focus on one 
random continuous initial parameter.


\noindent{\em Notation.}
For an interval $[r]=[\underline{r}, \overline{r}] \subset \mathbb{R}$ we denote the size of
the interval by $width([r]) = \overline{r} - \underline{r}$ and by                         
$mid([r]) = \frac{\overline{r} + \underline{r}}{2}$ the central point of the interval.

%\subsection{Interval Integration}
Our validated integration procedure employs the (1/3) Simpson rule:
\begin{equation}\label{eq:Simpson}
\begin{split}
K([I])= \int_{a}^{b} \, f(x)\, dx\, = \frac{width([I])}{6}(f(\underline{I}) + 4 f(mid([I])) + \\
f(\overline{I})) - \frac{width([I])^5}{2880}f^{(4)}(\xi)
\end{split}
\end{equation}
where $[I]=[a,b]$, $\xi \in [I]$ and $f^{(4)}$ is the fourth derivative of an integrable 
function $f$. For our applications the integrands are probability density functions, which
satisfy the required integrability and differentiability conditions.
Our aim is to compute an interval of arbitrary size $\epsilon \in (0, 1]\cap \mathbb{Q}$ 
that contains $K$.
 
\begin{definition}
An interval extension of function $f:X \rightarrow Y$ is an operator $[\cdot]$ such that:
\begin{equation*}
	\forall x \in [r] \subseteq X: f(x) \in [f]([r]) \subseteq Y
\end{equation*}
\end{definition} 
By applying interval arithmetics, one computes interval extensions of $f$ and $f^{(4)}$.
(Interval extensions can be computed using interval arithmetics libraries, \eg, FILIB++ \cite{filib}.)
The interval version of Simpson's rule can be obtained simply by replacing in (\ref{eq:Simpson})
the occurrences of $f$ and $f^{(4)}$ with their interval extension $[f]$, and by replacing $\xi$ 
with the entire interval $I$ \cite{ValidatedIntegration}:
\[
\begin{split}
K \in [K]([I]) = \frac{width([I])}{6}([f](\underline{I}) + 4 [f](mid([I])) + 
[f](\overline{I})) - \frac{width([I])^5}{2880}[f]^{(4)}([I]).
\end{split}
\]
Furthermore, by the definition of integral:
\begin{equation}\label{eq:K}
K \in \Sigma_{i = 1}^{n} [K]([x]_{i})
\end{equation}
where the collection of $[x]_{i}$'s is a partition of $[a, b]$. Note that we require a partition in a 
measure-theoretic sense, \ie, intersections have (Lebesgue) measure 0, since these have no effect on integration.

In order to guarantee $\epsilon$-completeness of the integration it is sufficient to partition 
$[a, b]$ into $n$ intervals $[x]_{i}$ such that for each $[x]_{i}$ we have
$width([I]([x]_{i})) < \epsilon \frac{width([x]_{i})}{b - a}$. Then, the exact value $K$ of the 
integral will belong to an interval (\ref{eq:K}) of width smaller than $\epsilon$. Pseudo-code 
for the procedure computing integral (\ref{eq:Simpson}) up to an arbitrary 
$\epsilon \in (0,1]\cap\mathbb{Q}$ is given in Algorithm \ref{alg:integration}. For our
purposes we will only make use of the interval partition $T$, which will enable us 
to compute the reachability probability, \ie, integral (\ref{eq:probabilistic-reachability_rv}),
with precision $\epsilon$.
\begin{proposition}\label{prop:IntegrationComplexity}
If $f \in P_{C^5[a,b]}$, then the complexity of Algorithm \ref{alg:integration} is {\em NP}.
\end{proposition}
\begin{algorithm}
\label{alg:integration}
input: function $f$, interval $[a, b]$, $\epsilon \in \mathbb{Q}^+$\;
output: $[I]$, partition $T$ of $[a, b]$ such that $\int_a^b f \in [I]$ 
	and $width([I]) \leq \epsilon$\;
$[I] = [0.0, 0.0]$\;
%$[K]([a,b]) = \frac{width([a, b])}{6}([f](a) + 4 [f](mid([a,b])) + [f](b)) - \frac{(width([a,b]))^5}{2880}[f^{(4)}]([a,b])$\;
$T,B = \emptyset$\;
\tcp{put initial partition on a stack}
$B.push(\{[a, b], [K]([a,b])\})$\;	
\While{$size(B) > 0$}{
	%$\{[x], [K]([x])\} = B.pop()$\;
	$\{[x], [y]\} = B.pop()$\;
	%\If{$width([K]([x])) > \epsilon \frac{width([x])}{b - a}$}{
	\If{$width([y]) > \epsilon \frac{width([x])}{b - a}$}{
		\tcp{split the interval in two}
		$B.push(\{[\underline{x}, mid([x])], [K]([\underline{x}, mid([x])])\})$\;
		$B.push(\{[mid([x]), \overline{x}], [K]([mid([x]), \overline{x}])\})$\;
	}\Else{
		\tcp{add sub-integral to the partial sum; save interval}
		$[I] = [I] + [K]([x])$\;
		%$T.push(\{[x], [K]([x])\})$\;
		$T.push([x])$\;
	}
}
\Return $T$, $[I]$\;
\caption{Validated Integration Procedure}
\end{algorithm}

\hide{
We recall that $P_{C^n[a,b]}$ denotes the class of polynomial-time (Type 2) computable functions
whose derivative $f^{(n)}$ exists and is continuous over $[a,b]$. Ko \cite[Section 6.2]{kobook}
showed that if $f$ is also analytic, then integration becomes P. However, such an algorithm
essentially uses truncated Taylor series over an arbitrary partition. Instead, 
Algorithm \ref{alg:integration} adaptively searches for a partition that guarantees the
required error bound $\epsilon$, while having a minimal number of intervals. In practice, this 
significantly benefits the performance of our whole implementation. Another advantage of 
Algorithm \ref{alg:integration} is that it does not require $f$ to be analytic.
}



results.tex
\vspace{-3ex}
\section{Experiments}\label{sec:Exp}
\vspace{-1ex}

We have implemented our algorithms in \preach; its source code 
and the models studied are on \texttt{https://github.com/dreal/probreach}. 
(The tool implementation is explained in \cite{ProbReach}.) 
The results below can be also accessed on 
\texttt{https://homepages.ncl.ac.uk/f.shmarov/probreach}. All experiments were
carried out on a multi-core Intel Xeon E5-2690 2.90GHz system running Linux Ubuntu 14.04LTS. 
The algorithms were also parallelised, and the results below feature the \preach\  CPU time of the
parallel version on 24 cores.

We have applied \preach\ to four hybrid models: a 2D-moving bouncing ball, human starvation, prostate
cancer therapy, and car collision scenario. The models feature a variety of highly nontrivial dynamics.
For example, the ODEs for the prostate cancer therapy model \cite{Liu2015} include exponential terms:
\[
\begin{split}
\frac{dx}{dt} = &\left(\frac{\alpha_{x}}{1+\mathrm{e}^{(k_1-z)k_2}}-\frac{\beta_{x}}{1+\mathrm{e}^{(z-k_3)k_4}}-m_1\left(1-\frac{z}{z_0}\right)-c_1\right)x+c_2 \\
\frac{dy}{dt} = &m_1\left(1-\frac{z}{z_0}\right)x + \left(\alpha_{y}\left(1-d_0\frac{z}{z_0}\right)-\beta_{y}\right)y \\
\frac{dz}{dt} = &-z \gamma - c_3
\end{split}
\]
More details on the models used in the experiments and the actual \preach\ model file for the
prostate cancer therapy can be found in Appendix \ref{apndx:Models}.

All experiments were validated using Monte Carlo probability estimation in MATLAB (the
reported CPU times are for one core).
In particular, we calculated confidence intervals using the 
%Chernoff-Hoeffding \cite{hoeffding} bound $N = \frac{\log{\frac{1}{1-c}}}{2\zeta^2}$, 
%where $\zeta$ is the interval half-width and $c$ is the coverage probability.
Chernoff-Hoeffding bound \cite{hoeffding}.
All results are given in the tables below, where the top half of each table contains the results 
obtained using our approach (\preach), while the bottom half reports the Monte Carlo results. 
For space reasons, the results of the bouncing ball are presented in the Appendix \ref{apndx:Models}.
Monte Carlo simulation of continuous nondeterminism (for \textbf{NPHA} models)
was implemented by first uniformly discretising the domain of the nondeterministic 
parameters. Then, for each (discretised) 
value of the parameters we built a confidence interval using the Chernoff-Hoeffding bound. 
Finally, the Monte Carlo interval reported in the tables below is the union of all such confidence 
intervals. Note that Monte Carlo intervals for \textbf{NPHA}s will be in general larger than $2\zeta$, and 
they will not be proper confidence intervals because of the nondeterministic parameters.
From the results we can see that our technique performs well even on highly nonlinear ODEs models such
as the prostate cancer treatment model, despite having unavoidably high complexity (see 
Theorem \ref{thm:complexity}). All Monte Carlo 
intervals cover the enclosures computed by \preach, thus confirming the correctness of our algorithms 
and their implementation. 

%A bigger variance for Gaussian parameters will not change the complexity significantly. The number of intervals will decrease: they will become larger and more evenly distributed over the domain (instead of having a large number of small intervals concentrated around the mean). We used small variances to investigate the performance of our approach on steep functions (such as a Gauss density with small variance), which are usually harder to integrate.
\vspace{-1ex}

\paragraph{\bf Table legend:}
{\small
$k$ = number of discrete transitions; $\epsilon$ = desired size of probability interval 
(\textbf{PHA}s only; lower box size limit for \textbf{NPHA});
$length$ = length of probability interval returned by \preach;
$\zeta,c$ = half-interval width and coverage probability for Chernoff bound;
$N$ = sample size from Chernoff bound;
$CPU$ = CPU time (sec).
}

\hide{
\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\centering
\caption{Bouncing ball model; see legend in Section \ref{sec:Exp}.}
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method & $k$ & $\epsilon$ & $length$ & & Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 
\hline 
\multirow{4}{4em}{{\bf Prob Reach}} & 0 & $10^{-9}$ & $5\cdot10^{-10}$ & &[8.21757e-05, 8.21762e-05] & 64 & 5\\ [0.5ex]
& 1 & $10^{-9}$ & $10^{-9}$ & &[0.1379483631, 0.1379483641] & 192 & 29\\ [0.5ex] 
%& 2 & $10^{-9}$ & $9.9\cdot10^{-10}$ & &[0.50868960502, 0.50868960601] & 927 & 164\\ [0.5ex] 
& 3 & $10^{-9}$ & $8\cdot10^{-10}$ & &[0.7387674005, 0.7387674013] & 3,806 & 563\\ [0.5ex] 
\hline \hline
Method & $k$ & $\zeta$ & $c$ & $P$ & Monte Carlo interval & $CPU_{seq}$ & $N$\\ [0.5ex] 
\hline 
\multirow{4}{4em}{Monte Carlo}& 0 & $5\cdot 10^{-6}$ & 0.99999 & 8.220032e-05 & [7.720032e-05, 8.720032e-05] & 16,455 & 230,258,509,300\\ [0.5ex]
& 1 & $5\cdot10^{-6}$ & 0.99999 & 0.1379449 & [0.1379399, 0.1379499] & 19,646 & 230,258,509,300 \\ [0.5ex] 
%& 2 & $5\cdot10^{-6}$ & 0.99999 & 0.5086939 & [0.5086889, 0.5086989] & 21,197 & 230,258,509,300 \\ [0.5ex] 
& 3 & $5\cdot10^{-6}$ & 0.99999 & 0.7387684 & [0.7387634, 0.7387734] & 20,975 & 230,258,509,300 \\ [0.5ex] 
\hline 
\end{tabular} 
\label{table:bouncing-ball} 
\end{adjustwidth}
\end{table} 
}

\hide{
\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\centering
\caption{Bouncing ball model; see legend in Section \ref{sec:Exp}.}
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method & Model type& $k$ & $\epsilon$ & $length$ & Probability interval & $CPU$ & \\ [0.5ex] 
\hline 
\multirow{2}{4em}{{\sffamily Prob Reach}} &\textbf{NPHA}& 0 & $10^{-3}$ & $2.38\cdot10^{-4}$ & [0.000013103, 0.000250681] & 223 & \\ [0.5ex]
&\textbf{NPHA}& 1 & $10^{-3}$ & $6.464\cdot10^{-2}$ &[0.0647381, 0.12937951] & 1,605 & \\ [0.5ex] 
\hline \hline
Method & Model type & $k$ & $\zeta$ & $c$ & Monte Carlo interval & $CPU$ & $N$\\ [0.5ex] 
\hline 
\multirow{2}{4em}{Monte Carlo}&\textbf{NPHA}& 0 & $5\cdot 10^{-3}$ & 0.99 & [0, 0.00520629] & 1,482 & 92,104\\ [0.5ex]
&\textbf{NPHA}& 1 & $5\cdot10^{-3}$ & 0.99 & [0.0585, 0.1367] & 1,485 & 92,104 \\ [0.5ex] 
\hline 
\end{tabular} 
\label{table:bouncing-ball} 
\end{adjustwidth}
\end{table} 
}

\hide{
\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\centering
\caption{Controlled bouncing ball model; see legend in Section \ref{sec:Exp}.}
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method&$k$ & $\epsilon$ & $length$ & &Probability interval & $CPU_{seq}$ & $CPU_{par}$\\ [0.5ex] 

\hline 
\multirow{2}{4em}{{\bf Prob Reach}}&2 & $10^{-2}$ & $8\cdot10^{-3}$ & &[0.199, 0.207] & 70 & 15\\ [0.5ex]
&2 & $10^{-9}$ & $3\cdot10^{-10}$ & &[0.2049030217, 0.204903022] & 8,332 & 1,156\\ [0.5ex] 
\hline\hline  
Method &$k$ & $\zeta$ & $c$ & $P$ & Monte Carlo interval & $CPU_{seq}$ & $N$\\ [0.5ex] 
\hline 
Monte Carlo&2 & $5\cdot10^{-3}$ & 0.99 & 0.2045948 & [0.1995948, 0.2095948] & 50,528 & 92,104\\ [0.5ex]
\hline 
\end{tabular} 
\label{table:controlled-bouncing-ball} 
\end{adjustwidth}
\end{table} 
}

\vspace{-3ex}

\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\caption{Starvation model; see legend in Section \ref{sec:Exp}.}
\centering
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method&Model type&$k$ & $\epsilon$ & $length$ & Probability interval & $CPU$ &\\ [0.5ex] 

\hline 
\multirow{2}{4em}{{\sffamily Prob Reach}}&\textbf{NPHA}&0 & $10^{-3}$ & $4.245\cdot10^{-3}$ &  [0.9219413, 0.92618671] & 1,152 & \\ [0.5ex]
&\textbf{PHA}&0 & $10^{-3}$ & $6.795\cdot10^{-4}$ &  [0.92455817, 0.92523768] & 23 & \\ [0.5ex]
\hline\hline  
Method &Model type&$k$ & $\zeta$ & $c$ &  Monte Carlo interval & $CPU$ & $N$\\ [0.5ex] 
\hline 
\multirow{2}{4em}{Monte Carlo}&\textbf{NPHA}& 0 & $5\cdot 10^{-3}$ & 0.99 &  [0.9179, 0.9311] & 12,433 & 92,104\\ [0.5ex]
&\textbf{PHA}& 0 & $5\cdot10^{-3}$ & 0.99 &  [0.9193355, 0.9293355] & 2,868 & 92,104 \\ [0.5ex] 
\hline

\end{tabular} 
\label{table:starvation} 
\end{adjustwidth}
\end{table} 

\vspace{-8ex}

\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\caption{Prostate cancer therapy model; see legend in Section \ref{sec:Exp}.}
\centering
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method&Model type&$k$ & $\epsilon$ & $length$ & Probability interval & $CPU$ &\\ [0.5ex] 

\hline 
\multirow{2}{4em}{{\sffamily Prob Reach}}&\textbf{PHA}&1 & $10^{-3}$ & $6.022\cdot10^{-4}$ & [0.47380981, 0.47441201] & 737 & \\ [0.5ex]
&\textbf{NPHA}&1 & $10^{-4}$ & $1.763\cdot10^{-3}$ & [0.4725522, 0.47431526] & 89,925 & \\ [0.5ex]
\hline\hline  
Method& Model type &$k$ & $\zeta$ & $c$ & Monte Carlo interval & $CPU$ & $N$\\ [0.5ex] 
\hline 
\multirow{2}{4em}{Monte Carlo}&\textbf{PHA}&1 & $1\cdot10^{-2}$ & 0.99 & [0.4648111, 0.4848111] & 5,700 & 23,026 \\ [0.5ex] 
&\textbf{NPHA}&1 & $1\cdot 10^{-2}$ & 0.99 & [0.4583, 0.4890] & 12,309 & 23,026\\ [0.5ex]
\hline

\end{tabular} 
\label{table:cancer} 
\end{adjustwidth}
\end{table} 

\vspace{-2ex}
\begin{table}[ht!] 
\begin{adjustwidth}{-1in}{-1in}
\caption{Car collision model; see legend in Section \ref{sec:Exp}.}
\centering
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method&Model type&$k$ & $\epsilon$ & $length$ & Probability interval & $CPU$ &\\ [0.5ex] 

\hline 
{\preach}&\textbf{PHA}&4 & $10^{-3}$ & $8.369\cdot10^{-4}$ & [0.5063922, 0.5072291] & 1,869 & \\ [0.5ex]
\hline\hline  
Method &Model type&$k$ & $\zeta$ & $c$ & Monte Carlo interval & $CPU$ & $N$\\ [0.5ex] 
\hline 
Monte Carlo&\textbf{PHA}&4 & $5\cdot 10^{-3}$ & 0.99 & [0.496629, 0.506629] & 32,201 & 92,104\\ [0.5ex]
\hline

\end{tabular} 
\label{table:cars} 
\end{adjustwidth}
\end{table} 

\hide{
\begin{table}[ht!] 
\caption{Insulin-glucose regulatory model; see legend in Section \ref{sec:Exp}. In this 
table, \textbf{ProbReach} $CPU_{seq}$ is the sum of all the actual CPU times of the parallel 
implementation.}
\begin{adjustwidth}{-1in}{-1in}
\centering
\vspace{-1ex}
\begin{tabular}{c c c c c c c r}
\hline
Method&$k$ & $\epsilon$ & $length$ && Probability interval & $CPU_{seq}$ & $CPU_{par}$\\
\hline 
\multirow{3}{4em}{{\bf Prob Reach}}&0 & $10^{-2}$ & $5.328\cdot10^{-3}$ & &[0.994589, 0.999917] & 2,805,634 & 165,404\\ [0.5ex]
&1 & $10^{-3}$ & $8.1\cdot10^{-4}$ & &[0.999107, 0.999917] & 3,326,581 & 443,910\\ [0.5ex]
&1 & $10^{-4}$ & $5.5\cdot10^{-5}$ & &[0.999657, 0.999712] & 3,498,765 & 490,257\\ [0.5ex] 
\hline\hline 
Method &$k$ & $\zeta$ & $c$ & $P$ & Monte Carlo interval & $CPU_{seq}$ & $N$\\ [0.5ex] 
\hline
\multirow{2}{4em}{Monte Carlo}&1 & $5\cdot10^{-3}$ & 0.99 & 0.997266555 &[0.9945331, 1] & 58,069 & 92,104\\ [0.5ex]
&1 & $2.5\cdot10^{-3}$ & 0.99 & 0.99853 &[0.99706, 1] & 219,623 & 368,416\\ [0.5ex]
%\multirow{2}{4em}{Monte Carlo}&1 & $5\cdot10^{-3}$ & 0.99 & [0.9945331, 1.0045331] & 58,069 & -\\ [0.5ex]
%&1 & $2.5\cdot10^{-3}$ & 0.99 & [0.997052137, 1.002052137] & 404,495 & -\\ [0.5ex]
\hline
\end{tabular} 
\label{table:insulin-glucose} 
\end{adjustwidth}
\end{table} 
}

%We haven't included the IG model, but perhaps a similar thing holds for another model?
\hide{
Second, except for the insulin-glucose model 
(Table \ref{table:insulin-glucose}), ProbReach is orders of magnitude faster than Monte Carlo, 
while providing shorter and fully verified probability intervals. The insulin-glucose
model is a demanding benchmark, and while ProbReach clearly loses for $\epsilon = 10^{-2}$, we can
see from Table \ref{table:insulin-glucose} that obtaining a Monte Carlo confidence interval of size
$5.5\times 10^{-5}$ would require about $58,069 \times 10,000 = 580,690,000$ seconds. (With the
Chernoff bound a 10-fold reduction of $\zeta$ increases the sample size 100-fold.) ProbReach
computed such an interval in about 3,500,000 seconds.
}


main.bbl
\begin{thebibliography}{10}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }

\bibitem{Abate2010624}
Abate, A., Katoen, J.P., Lygeros, J., Prandini, M.: Approximate model checking
  of stochastic hybrid systems. European Journal of Control  16(6),  624 -- 641
  (2010)

\bibitem{DBLP:conf/hybrid/AlurCHH92}
Alur, R., Courcoubetis, C., Henzinger, T.A., Ho, P.H.: Hybrid automata: An
  algorithmic approach to the specification and verification of hybrid systems.
  In: Hybrid Systems. LNCS, vol. 736, pp. 209--229 (1992)

\bibitem{DBLP:conf/icalp/AlurD90}
Alur, R., Dill, D.L.: Automata for modeling real-time systems. In: ICALP. LNCS,
  vol. 443, pp. 322--335 (1990)

\bibitem{BrihayeDGORW11}
Brihaye, T., Doyen, L., Geeraerts, G., Ouaknine, J., Raskin, J., Worrell, J.:
  On reachability for hybrid automata over bounded time. In: ICALP. LNCS, vol.
  6756, pp. 416--427 (2011)

\bibitem{UPPAALSMCtut}
David, A., Larsen, K., Legay, A., Miku\v{c}ionis, M., Poulsen, D.B.: Uppaal
  {SMC} tutorial. International Journal on Software Tools for Technology
  Transfer (STTT)  (2015), {\em to appear.}

\bibitem{SMCnd}
Ellen, C., Gerwinn, S., Fr{\"{a}}nzle, M.: Statistical model checking for
  stochastic hybrid systems involving nondeterminism over continuous domains.
  International Journal on Software Tools for Technology Transfer (STTT)
  (2014), {\em to appear.}

\bibitem{Enszer10verifiedsolution}
Enszer, J.A., Stadtherr, M.A.: Verified solution and propagation of uncertainty
  in physiological models. Reliable Computing  15,  168--178 (2010)

\bibitem{FraenzleHHWZ11}
Fr{\"a}nzle, M., Hahn, E.M., Hermanns, H., Wolovick, N., Zhang, L.:
  Measurability and safety verification for stochastic hybrid systems. In:
  HSCC. pp. 43--52 (2011)

\bibitem{Sisat}
Fr{\"{a}}nzle, M., Teige, T., Eggers, A.: Engineering constraint solvers for
  automatic analysis of probabilistic hybrid automata. J. Log. Algebr. Program.
   79(7),  436--466 (2010)

\bibitem{ValidatedIntegration}
Galdino, S.: Interval integration revisited. Open Journal of Applied Sciences
  2(4B),  108--111 (2012)

\bibitem{DBLP:conf/lics/GaoAC12}
Gao, S., Avigad, J., Clarke, E.M.: Delta-decidability over the reals. In: LICS.
  pp. 305--314 (2012)

\bibitem{DBLP:journals/corr/GaoKCC14}
Gao, S., Kong, S., Chen, W., Clarke, E.M.: Delta-complete analysis for bounded
  reachability of hybrid systems. CoRR  arXiv:1404.7171 (2014)

\bibitem{DBLP:conf/cade/GaoKC13}
Gao, S., Kong, S., Clarke, E.M.: {dReal}: An {SMT} solver for nonlinear
  theories over the reals. In: CADE-24. LNCS, vol. 7898, pp. 208--214 (2013)

\bibitem{hoeffding}
Hoeffding, W.: Probability inequalities for sums of bounded random variables.
  J. Amer. Statist. Assoc.  58(301),  13--30 (1963)

\bibitem{Ideta2008}
Ideta, A.M., Tanaka, G., Takeuchi, T., Aihara, K.: A mathematical model of
  intermittent androgen suppression for prostate cancer. Journal of Nonlinear
  Science  18(6),  593--614 (2008)

\bibitem{kobook}
Ko, K.I.: Complexity Theory of Real Functions. Birkh\"{a}user (1991)

\bibitem{KNP11}
Kwiatkowska, M., Norman, G., Parker, D.: {PRISM} 4.0: Verification of
  probabilistic real-time systems. In: CAV. LNCS, vol. 6806, pp. 585--591
  (2011)

\bibitem{UPPAAL}
Larsen, K.G., Pettersson, P., Yi, W.: Uppaal in a nutshell. International
  Journal on Software Tools for Technology Transfer (STTT)  1,  134--152 (1997)

\bibitem{filib}
Lerch, M., Tischler, G., Gudenberg, J.W.V., Hofschuster, W.e., Kr\"{a}mer, W.:
  {FILIB++}, a fast interval library supporting containment computations. ACM
  Trans. Math. Softw.  32(2),  299--324 (2006)

\bibitem{Liu2015}
Liu, B., Kong, S., Gao, S., Zuliani, P., Clarke, E.M.: Towards personalized
  cancer therapy using delta-reachability analysis. In: HSCC. pp. 227--232. ACM
  (2015)

\bibitem{Novak1992}
Novak, E., Wo\'{z}niakowski, H.: Relaxed verification for continuous problems.
  Journal of Complexity  8(2),  124 -- 152 (1992)

\bibitem{Petras:VerifiedNI}
Petras, K.: Principles of verified numerical integration. Journal of
  Computational and Applied Mathematics  199(2),  317 -- 328 (2007)

\bibitem{RamponiCSL10}
Ramponi, F., Chatterjee, D., Summers, S., Lygeros, J.: On the connections
  between {PCTL} and dynamic programming. In: HSCC. pp. 253--262. {ACM} (2010)

\bibitem{Richardson68}
Richardson, D.: Some undecidable problems involving elementary functions of a
  real variable. J. Symb. Log.  33(4),  514--520 (1968)

\bibitem{ProbReach}
Shmarov, F., Zuliani, P.: {ProbReach}: Verified probabilistic
  $\delta$-reachability for stochastic hybrid systems. In: HSCC. pp. 134--139.
  ACM (2015)

\bibitem{Starvation}
Song, B., Thomas, D.: Dynamics of starvation in humans. Journal of Mathematical
  Biology  54(1),  27--43 (2007)

\bibitem{faust2}
Soudjani, S.E.Z., Gevaerts, C., Abate, A.: {FAUST}\({}^{\mbox{2}}\): Formal
  abstractions of uncountable-state stochastic processes. In: TACAS (2015),
  {\em to appear.}

\bibitem{Tiwari02}
Tiwari, A.: Formal semantics and analysis methods for {S}imulink {S}tateflow
  models. Tech. rep., SRI International (2002),
  \url{http://www.csl.sri.com/users/tiwari/html/stateflow.html}

\bibitem{TkachevA13}
Tkachev, I., Abate, A.: Formula-free finite abstractions for linear temporal
  verification of stochastic hybrid systems. In: HSCC. pp. 283--292. {ACM}
  (2013)

\bibitem{ZhangSRHH10}
Zhang, L., She, Z., Ratschan, S., Hermanns, H., Hahn, E.M.: Safety verification
  for probabilistic hybrid systems. In: CAV. LNCS, vol. 6174, pp. 196--211
  (2010)

\end{thebibliography}
