text:
\section{Acknowledgments}
\label{sec:acknowledgments}

The authors gratefully acknowledge partial financial supports from the
research training group CONVEY funded by the German Research Foundation under grant GRK 2428.

\bibliographystyle{ieeetr}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we presented a novel generic framework
to automatically tune all algorithm parameters{}
which is a major problem of present reachability algorithms.
The presented algorithm enables a fully-automated computation of the reachable set
whose error is below a user-defined error.
Previous work only considered the tuning of time parameters.
An example implementation has shown to outperform manually-tuned algorithm parameters{} on benchmarks
and provides better results than genetic algorithms searching for algorithm parameters{}
on randomly-generated systems of varying dimensions.
The extension of the presented framework to nonlinear systems
will be considered in the future.

\section{Implementation}
\label{sec:implementation}

In this section, we present an implementation of our framework
and validate our assumptions using a concrete reachability algorithm and a concrete set representation.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Computation of Reachable Sets}
\label{ssec:setprop}

We over-approximate the exponential matrix~$e^{At}$
by a finite number of Taylor terms~$\eta \in \naturalsof{+}$
and an interval matrix~$\mathbf{E}{}$ enclosing the remainder \cite[Prop. 2]{Althoff2011a}.
\begin{align}
	e^{A \Delta t} &\in \sum_{k=1}^{\eta} \frac{(A \Delta t)^k}{k!} \oplus \Eof{\Delta t, \eta} \; ,
		\label{eq:eAt} \\
	\Eof{\Delta t, \eta} &= [-\Eabsof{\Delta t, \eta}, \Eabsof{\Delta t, \eta}]\label{eq:E}  \\
	\text{with} \quad \Eabsof{\Delta t, \eta}
		&= \Big| \sum_{k=\eta + 1}^\infty \frac{1}{k!} \big(|A| \Delta t\big)^k \Big| \; . \label{eq:Eabs}
\end{align}
To cover all trajectories in a time interval spanned by~$\Delta t$,
we introduce the terms $\mathbf{F}_x{}$ \cite[Sec. 4]{Althoff2011a}
and $\mathbf{F}_u{}$ \cite[Sec. 3.2.2]{Althoff2010a}:
\begin{align}
	\Fof{\Delta t, \eta} &= \bigoplus_{k=2}^\eta [(k^{\frac{-k}{k-1}}
		- k^{\frac{-1}{k-1}}) \Delta t^k,0] \frac{A^k}{k!}
		\oplus \Eof{\Delta t, \eta} \; , \label{eq:F} \\
	\begin{split} \label{eq:inputCorr}
		\inputCorrof{\Delta t, \eta} &= \bigoplus_{k=2}^{\eta+1}
			[(k^{\frac{-k}{k-1}} - k^{\frac{-1}{k-1}}) \Delta t^k,0] \frac{A^{k-1}}{k!} \\
			&\qquad \oplus \Eof{\Delta t, \eta} \, \Delta t \; .
	\end{split}
\end{align}

The reachable sets for the homogeneous and inhomogeneous solution,
generally defined in~\eqref{eq:HRtieqplus} and~\eqref{eq:PRtieqplus},
are computed by~\cite[Sect. 3.2.1-3.2.2]{Althoff2010a}
\begin{align}
	\HRti{[0, \Delta t]} &=
		\underbrace{\CH{\mathcal{X}^0{}, e^{A \Delta t} \mathcal{X}^0{}}}_{= \, \HReqtiunderbrace{[0,\Delta t]}} \nonumber \\
		&\qquad \oplus \underbrace{\Fof{\Delta t, \eta} \, \mathcal{X}^0{}
		\oplus \inputCorrof{\Delta t, \eta} \, c_{u}{}
			_{= \, \HRplustiunderbrace{[0,\Delta t]}} \;, \label{eq:HRti_cora} \\
	\PRti{[0, \Delta t]} &= \underbrace{\sum_{k=0}^{\eta}
		\Big( \frac{A^k \Delta t^{k+1}}{(k+1)!} \Big) \, \mathcal{U}{}}_{= \, \PReqtiunderbrace{[0,\Delta t]}}
		\oplus \underbrace{\Eof{\Delta t, \eta} \, \Delta t \, \mathcal{U}{}}_{= \, \PRplustiunderbrace{[0,\Delta t]}} ,
			\label{eq:PRti_cora}
\end{align}
with~$c_{u}{}$ being the center of the input set~$\mathcal{U}{}$.
We include the term~$\inputCorrof{\Delta t, \eta} \, c_{u}{}$ in the homogeneous solution
as it only covers the current time interval and therefore does not accumulate over time.
In the next section, we will verify the applicability of
\eqref{eq:HRti_cora} and \eqref{eq:PRti_cora} for our adaptive framework.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Verification of Assumptions}
\label{ssec:verass}

We now want to verify \assref{ass:bloatHRPR} and \propref{prop:bloatPRtiadm}.
To obtain the error $\bloatHRti{[t_i, t_{i+1}]}$
of the homogeneous solution $\HRti{[t_i, t_{i+1}]}$,
we multiply \eqref{eq:HRti_cora} by $e^{A t_i}$ as required in \eqref{eq:HRPR_start}
and insert the result in \eqref{eq:bloatHRti}, which yields
\begin{align}
	\begin{split} \label{eq:bloatHRti_cora}
	\bloatHRti{[t_i, t_{i+1}]} &=
		\errOp{ e^{A t_i} \, \Fof{\Delta t_i, \eta_i} \, \mathcal{X}^0{} \\
			&\qquad \qquad \oplus \, e^{A t_i} \, \inputCorrof{\Delta t_i, \eta_i} \, c_{u}{} } \; .
	\end{split}
\end{align}
\begin{proposition} \label{prop:bloatHRti_cora}
The error $\bloatHRti{[t_i,t_{i+1}]}$ in \eqref{eq:bloatHRti_cora}
converges to 0 for $\Delta t_i \rightarrow 0$
and therefore satisfies \assref{ass:bloatHRPR}.
\end{proposition}
\begin{proof}
Since the operation $\errOp{\mathcal{S}{}}$ returns the enclosing radius
of the interval over-approximation of a set $\mathcal{S}{}$,
it suffices to show that the volume of $\mathcal{S}{}$
converges to 0 for $\Delta t_i \rightarrow 0$.
This is the case if the interval matrices $\mathbf{F}_x{}$ and $\mathbf{F}_u{}$
converge to $[\zeroMat{n}, \zeroMat{n}]$:
Using \eqref{eq:inputCorr}, we yield
$\lim_{\Delta t \rightarrow 0} \inputCorrof{\Delta t, \eta} = [\zeroMat{n}, \zeroMat{n}]$.
For $\mathbf{F}_x{}$, we have that
$\lim_{\Delta t \rightarrow 0} \Eof{\Delta t, \eta} \overset{\eqref{eq:E}}{=} [\zeroMat{n},\zeroMat{n}]$ since
$\lim_{\Delta t \rightarrow 0} \Eabsof{0, \eta} \overset{\eqref{eq:Eabs}}{=} \zeroMat{n}$.
By plugging this in \eqref{eq:F}, we immediately see that
$\Fof{\Delta t, \eta} = [\zeroMat{n}, \zeroMat{n}]$ for $\Delta t \rightarrow 0$.
\end{proof}

To obtain the error $\bloatPRti{[t_i, t_{i+1}]}$
of the inhomogeneous solution $\PRti{[t_i, t_{i+1}]}$,
we multiply \eqref{eq:PRti_cora} by $e^{A t_i}$ as required in \eqref{eq:HRPR_start}
and insert the result in \eqref{eq:bloatPRti}, which gives us
\begin{equation} \label{eq:bloatPRti_cora}
	\bloatPRti{[t_i, t_{i+1}]} =
		\errOp{ e^{A t_i} \, \Eof{\Delta t_i, \eta_i} \, \Delta t_i \, \mathcal{U}{} } \; .
\end{equation}
\begin{proposition} \label{prop:bloatPRti_cora}
The error $\bloatPRti{[t_i,t_{i+1}]}$ in \eqref{eq:bloatPRti_cora}
converges to 0 for $\Delta t_i \rightarrow 0$
and therefore satisfies \assref{ass:bloatHRPR}.
\end{proposition}
\begin{proof}
Equivalently to the proof of \propref{prop:bloatHRti_cora},
we show that the volume of the resulting set converges to 0 for ${\Delta t_i \rightarrow 0}$.
Since $\Delta t_i$ appears as a multiplicative factor, it holds that
$\lim_{\Delta t_i \rightarrow 0} \, e^{A t_i} \, \Eof{\Delta t_i, \eta_i} \, \Delta t_i \, \mathcal{U}{}
= [\zeroMat{n}, \zeroMat{n}]$.
$~$ \hfill $~$
\end{proof}
We introduce the following lemma for the subsequent derivations:
\begin{lemma} \label{lmm:E}
The size of $\Eof{\Delta t, \eta}$ decreases superlinearly with respect to $\Delta t$:
\begin{equation*}
	\forall \varphi \in (0,1): \Eof{\varphi \Delta t, \eta} \subseteq \varphi \, \Eof{\Delta t, \eta} \; .
\end{equation*}
\end{lemma}

\noindent \textit{Proof.}
Using $\Eof{\Delta t, \eta} \overset{\eqref{eq:E}}{=} [-\Eabsof{\Delta t, \eta},\Eabsof{\Delta t, \eta}]$,
it is sufficient to show that each entry of $\Eabsof{\Delta t, \eta}$ decreases superlinearly
with respect to $\Delta t$:
\begin{align*}
	\hspace{20pt}
	&\forall \varphi \in (0,1): \Eabsof{\varphi \Delta t, \eta} \overset{\eqref{eq:Eabs}}{=}
	\Big| \sum_{k = \eta + 1}^{\infty} \frac{1}{k!} \big(|A| \varphi \Delta t\big)^k \Big| \\
	&\qquad \leq \varphi \, \Big| \sum_{k = \eta + 1}^{\infty} \frac{1}{k!} \big(|A| \Delta t\big)^k \Big|
	= \varphi \, \Eabsof{\Delta t, \eta} \; .
	\hspace{18pt} \square
\end{align*}

\begin{theorem}
The error $\bloatPRti{[t_i,t_{i+1}]}$ in \eqref{eq:bloatPRti_cora}
decreases superlinearly according to \defref{def:superlinear}:
\begin{equation*}
	\forall \varphi \in (0,1): \bloatPRti{[t_i, t_i + \varphi \Delta t_i]}
		\leq \varphi \, \bloatPRti{[t_i, t_i + \Delta t_i]} \; .
\end{equation*}
\end{theorem}

\textit{Proof:}
Since $\errOp{\cdot}$ is defined by the enclosing radius~$r$
according to \propref{prop:bloat},
it suffices to show that the enclosing radius decreases
superlinearly with respect to $\Delta t_i$:
\begin{equation*}
	\forall \varphi \in (0,1): r \big( \PRplusti{[t_i, t_i + \varphi \Delta t_i]} \big)
	\leq \varphi \, r \big( \PRplusti{[t_i, t_i + \Delta t_i]} \big) .
\end{equation*}
This condition is satisfied if
\begin{equation*}
	\forall \varphi \in (0,1): \PRplusti{[t_i, t_i + \varphi \Delta t_i]}
	\subseteq \varphi \, \PRplusti{[t_i, t_i + \Delta t_i]} \; .
\end{equation*}
Using the definition of $\PRplusti{[0, \Delta t]}$ in \eqref{eq:PRti_cora}
and \lmmref{lmm:E} it holds that
\begin{align*}
	&\PRplusti{[t_i, t_i + \varphi \Delta t_i]}
		\overset{\eqref{eq:PRti_cora}}{=} ( \varphi \, \Delta t_i ) \, e^{A t_i} \,
		\underbrace{\Eof{\varphi \Delta t_i, \eta_i}}_%
			{\overset{\text{\lmmref{lmm:E}}}{\subseteq} \,
			\varphi \, \Eof{\Delta t_i, \eta_i}} \, \mathcal{U}{} \subseteq \\
		&\varphi \, \big( \underbrace{\varphi \Delta t_i \, e^{A t_i} \, \Eof{\Delta t_i, \eta_i} \, \mathcal{U}{}}_%
			{\overset{\eqref{eq:PRti_cora}}{=} \,
			\varphi \, \PRplusti{[t_i, t_i + \Delta t_i]}} \big)
		\overset{\varphi \in (0,1)}{\subseteq} \varphi \, \PRplusti{[t_i, t_i + \Delta t_i]}. \;
	\square
\end{align*}

The cut-off value $\eta_{\text{max}}{}$ can be automatically obtained according to \cite{Bickart1968}
to truncate the power series in \eqref{eq:eAt}.
If $\eta_{\text{max}}{}$ does not satisfy the current error bounds,
we proceed to smaller values for $\Delta t$
which are guaranteed to eventually satisfy the error by \thmref{thm:alg:full}.

Following the above derivations, the presented implementation satisfies \thmref{thm:alg:full}.
Hence, it can be used to perform reachability analysis for LTI systems
while remaining below a user-defined over-approximation error.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Set Representation}
\label{ssec:setrepr}

It remains to choose a set representation.
The work in~\cite{Althoff2016c} shows that zonotopes are optimal
and that one should add support functions if the initial set is not a zonotope.
Since we only use zonotopes as initial sets, we only use them:

\begin{definition} \label{def:zonotopes} (Zonotopes)
Given a center ${c \in \mathbb{R}{n}}$
and an arbitrary number ${\gamma{} \in \mathbb{N}{}}$
of generator vectors ${g^{(1)}, ..., g^{(\gamma{}{})} \in \mathbb{R}{n}}$,
a zonotope is defined as \cite[Def. 1]{Girard2005}
\begin{equation*}
	\mathcal{Z} = \Big\{ x \in \mathbb{R}{n} \, \Big| \, x
		= c + \sum_{i=1}^{\gamma{}} \beta_i \cdot g^{(i)}, \, -1 \leq \beta_i \leq 1 \Big\} \; .
\end{equation*}
The order of a zonotope is $\rho{} = \frac{\gamma{}}{n}$. \hfill $\square$
\end{definition}

Following \defref{def:zonotopes}, we have that $\Phi_{\text{set}}{} = \rho{}$.
We iteratively decrease the zonotope order by decrementing
the total number of generators describing $\PRti{[0,t_{i+1}]}$
as shown in~\cite{Kopetzki2017},
which also provides us with $\bloatSetti{[t_i,t_{i+1}]}$.
Therefore, we define the operator
$\decrOp{\Phi_{\text{set}}{}}: n\gamma{} \gets n(\gamma{} - 1)$.
In the next section, we apply the presented implementation.

\section{Introduction}
\label{sec:Introduction}

Reachability analysis is one of the main techniques
to formally verify the correctness of mixed discrete/continuous systems:
It computes the set of reachable states of a system
for a set of uncertain initial states as well as uncertain inputs.
If the reachable set does not intersect the unsafe regions
defined by a given safety property, that property is satisfied.

Exact reachable sets can only be computed
for a limited class of systems \cite{LafPap2001}.
In the remaining cases, one calculates over-approximations
whose tightness is strongly influenced by algorithm parameters{};
unlike model parameters, these parameters have no relation to the considered model.

Due to large over-approximations resulting from poor algorithm parameters{},
reachability analysis may fail to prove a safety property
even though the property is satisfied by the exact reachable set.
We address this problem by proposing a novel generic framework
to automatically tune all algorithm parameters{} for reachability analysis
of linear continuous time-invariant systems during runtime
respecting a user-defined error bound.
Our framework can be used for different reachability algorithms
and for different set representations.

\paragraph{State of the Art}

Reachability algorithms for linear continuous systems
are mainly based on the propagation of reachable sets
for multiple time steps until reaching a fixed point or a user-defined time horizon.
Propagation-based techniques have been extensively investigated
\cite{Althoff2007c,Che2015,Gurung2019,Girard2006,BFFPSV2018,FreLeG2011}
using tools such as
\textit{CORA} \cite{Althoff2015a}, \textit{Flow*} \cite{Che2013},
\textit{XSpeed} \cite{Ray2015b}, \textit{SpaceEx} \cite{FreLeG2011},
and \textit{JuliaReach} \cite{BogFor2019}.
A further method is to compute reachable sets using simulations
\cite{Bak2017b,Duggirala2016}, which is implemented
in the tool \textit{HyLAA} \cite{Bak2017}.
The used set representation is another distinctive feature
besides the method to compute reachable sets.
Many different set representations have been researched,
including polytopes \cite{FreLeG2011},
zonotopes \cite{Girard2006}, ellipsoids \cite{Kurzhanski2000},
griddy polyhedra \cite{AsaDan2000}, star sets \cite{Bak2017b},
support functions \cite{Girard2008b},
and constrained zonotopes \cite{Scott2016}.

Previous approaches for automated parameter tuning
mainly focused on the time step size:
Numerical ODE solvers often compute different solutions in parallel
and decrease the time step size if the difference
between solutions exceeds a certain threshold
\cite{Lapidus1971,Ascher1994}.
Furthermore, several automated time step adaptation strategies
have been developed for \textit{guaranteed integration} methods
that enclose only a single trajectory rather than a set of trajectories
\cite{Kerbl1991,Rufeger1993,Nedialkov2000}.
For reachability analysis, only a few methods automatically set the time step size.
The approach in \cite{FreLeG2011} chooses the time step size
automatically in order to keep the error in a user-defined direction
below a user-defined bound.
Furthermore, the approach in \cite{PraVis2011}
automatically chooses suitable time step sizes
so that the Hausdorff distance between the exact reachable set
and the computed over-approximation stays below a certain threshold.
So far, there is no approach that considers all algorithm parameters{}.

\paragraph{Contributions}

We introduce a novel generic framework that automatically tunes
all algorithm parameters{} of a reachability algorithm
so that the over-approximation error stays below a user-defined threshold.
This generalized framework can be applied to
many different reachability algorithms
and we show that our self-parametrization approach always converges
if the reachability algorithm satisfies certain requirements.

\medskip
After introducing some preliminaries in \sectref{sec:prelim},
we present our generic framework in \sectref{sec:selfparam}.
In \sectref{sec:implementation}, we demonstrate the implementation
of our framework for a specific reachability algorithm.
Finally, the evaluation of our implementation on several numerical examples
in \sectref{sec:numex} demonstrates
the fully automated computation
of tight over-approximations of reachable sets
with only little computational overhead.
\section{Numerical Examples}
\label{sec:numex}

We have implemented our approach in \textit{MATLAB}.
To extensively test our approach, we perform several investigations:
First, we measure the computational overhead caused by our parameter tuning.
Next, \algref{alg:full} is evaluated on benchmark systems and compared to manually-tuned algorithm parameters{}.
Finally, we compare \algref{alg:full} with a genetic algorithm searching for algorithm parameters{}.
Since the genetic algorithm requires a lot of memory,
we used an Intel Xeon Gold 6136 3.00GHz processor and 768GB of DDR4 2666/3273MHz memory,
all other computations have been performed using an Intel i3 processor with 8GB memory.
For the genetic algorithm and the overhead measurement,
we used the subsequently introduced randomly-generated systems
because the results might vary depending on the investigated system.

\paragraph{Random Generation of Systems}
\label{par:randsys}

We randomly picked complex conjugate pairs of eigenvalues
from a uniform distribution over the range of $[-1,1]$ for the real part
and $[-\mathrm{i},\mathrm{i}]$ for the imaginary part
without loss of generality, since the characteristics of the solution
only depends on the ratio of real and imaginary values.
The state-space form for these eigenvalues is computed
and subsequently rotated by a random orthogonal matrix
of increasing sparsity for higher dimensions.
Furthermore, the initial set $\mathcal{X}^0{}$ and the input set $\mathcal{U}{}$
are given by hypercubes centered at $(10,...,10)^T$ with edge length 0.5,
and $(1,...,1)^T$ with edge length 0.1, respectively.
Lastly, we set $\varepsilon_{\text{max}}{} = 0.05$ and $T=3$s.

\begin{figure}[h]
\begin{center}
	\setlength{\belowcaptionskip}{-15pt}
	\includegraphics[width = 0.49 \textwidth]{Figures/Overhead_50sys_final.pdf}
	\caption{Time consumption by parameter tuning
		in relation to total execution time using
		50 randomly-generated systems per dimension.}
	\label{fig:overhead}
\end{center}
\end{figure}

\paragraph{Computational Overhead}
\label{par:compoverhead}

First, we investigate the computational overhead caused
by the continuous adaptation of all algorithm parameters{}.
For this purpose, we measured the time \algref{alg:full}
spends on the adaptation and the set propagation
over 50 randomly-generated systems of dimensions 5 to 75.
\figref{fig:overhead} shows the fraction of the total time
spent on the parameter tuning:
The overhead remains manageable even for high-dimensional systems,
settling at about the same time as used for the set propagation.
This is achieved by the computational efficiency of \propref{prop:bloat}.
The overhead is more than compensated by the adaptive adjustment of the algorithm parameters{}
as discussed at the end of this section.
Also, the common practice of trial and error requires tuning times
that exceed computation times by several factors.

\paragraph{ARCH benchmarks}
\label{par:ARCH}

Next, we have applied our algorithm to the
48-dimensional Building (BLD) benchmark
and the 273-dimensional International Space Station (ISS) benchmark,
both from the ARCH competition \cite{ARCH19}.
These systems are manually tuned for the competition
by executing many runs to optimize the algorithm parameters
with respect to the computation time
while still satisfying all specifications.
For a fair comparison, we set $\varepsilon_{\text{max}}{}$
to the highest possible value that still verifies the given specification.
We compare our results to the tool \textit{CORA} \cite{Althoff2015a}
since it is also implemented in \textit{MATLAB}.

\begin{table*}[t]
\begin{center}
	\caption{Results of ARCH benchmarks for \algref{alg:full} and \textit{CORA}.}
	\label{tab:ARCH}
	\begin{tabular}{l
					c @{\hspace{7pt}}
					c @{\hspace{7pt}}
					c @{\hspace{7pt}}
					c @{\hspace{7pt}}
					c
					c @{\hspace{7pt}}
					c @{\hspace{7pt}}
					c}
		\toprule
		\multirow{2}{2.1cm}{\textbf{Benchmark}} & \multicolumn{4}{c}{\textbf{\algref{alg:full}}} & &
			\multicolumn{3}{c}{\textbf{CORA}} \\ \cmidrule{2-5} \cmidrule{7-9}
		 & $\varepsilon_{\text{max}}{}$ & Time & Steps & $[\Delta t_{\text{min}}, \Delta t_{\text{max}}]$ &
		 & Time & Steps & $[\Delta t_{\text{min}}, \Delta t_{\text{max}}]$ \\ \cmidrule{1-9}
		BLDC01 & $2 \cdot 10^{-3}$ & 4.0s & 839 & $[0.0081,0.0448]$ & & 5.5s & 2400 & $[0.0020,0.0100]$ \\
		BLDF01 & $6 \cdot 10^{-3}$ & 5.6s & 818 & $[0.0081,0.0597]$ & & 6.0s & 2400 & $[0.0020,0.0100]$ \\
		ISSC01 & $5.6$ & 21.7s & 189 & $[0.0704,0.1371]$ & & 22.8s & 1000 & $[0.0200,0.0200]$ \\
		ISSF01 & $2 \cdot 10^{-3}$ & 295s & 1216 & $[0.0059,0.0395]$ & & 849s & 2000 & $[0.0100,0.0100]$ \\
		\bottomrule
	\end{tabular}
\end{center}
\vspace{-20pt}
\end{table*}

\tabref{tab:ARCH} shows a comparison between \algref{alg:full} and \textit{CORA}
in terms of the computation time, the number of steps,
and the minimum and maximum values for~$\Delta t$.
\algref{alg:full} adapts the values of the algorithm parameters{} depending on the current system behavior
by enlarging the time step size in regions where this only moderately increases the over-approximation.
The resulting range can be observed by the large differences
between $\Delta t_\text{min}$ and $\Delta t_\text{max}$.
Consequently, the total number of steps decreases making less computations
necessary than in the case of fixed algorithm parameters as used by \textit{CORA}
and shown in \tabref{tab:ARCH}.
The range for $\Delta t$ in both building benchmarks by \textit{CORA}
is explained by the switching between two manually-tuned time steps.

\begin{figure}[h]
\begin{center}
	\setlength{\belowcaptionskip}{-25pt}
	\includegraphics[width = 0.48 \textwidth]{Figures/ISSF01_diffeps_final.pdf}
	\caption{Benchmark ISSF01 with
		the specifications ISS01, ISU01,
		and the reachable sets $\Rti{[0,T]}$ in
		dark gray ($\varepsilon_{\text{max}}{} = 20 \cdot 10^{-3}$),
		light gray ($\varepsilon_{\text{max}}{} = 10 \cdot 10^{-3}$),
		and ivory ($\varepsilon_{\text{max}}{} = 2 \cdot 10^{-3}$).}
	\label{fig:diffbloatmax}
\end{center}
\end{figure}

\begin{figure*}
    \centering
    \setlength{\belowcaptionskip}{-15pt}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{Figures/ISSF01_tVec_final.pdf}
    \end{subfigure}
    ~

    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{Figures/ISSF01_etas_final.pdf}

    \end{subfigure}
    \begin{subfigure}[b]{0.3\textwidth}
        \includegraphics[width=\textwidth]{Figures/ISSF01_zonos_final.pdf}
    \end{subfigure}
    \captionsetup{
    	width=\textwidth,
    	labelfont=bf,
    	font=small,

	}
	\caption{Benchmark ISSF01: $\Delta t, \eta, \rho{}$ for different $\varepsilon_{\text{max}}{}$:
		black ($\varepsilon_{\text{max}}{}=20\cdot 10^{-3}$),
		gray ($\varepsilon_{\text{max}}{}=10\cdot 10^{-3}$),
		and blue ($\varepsilon_{\text{max}}{}=2\cdot 10^{-3}$).}
	\label{fig:ISSFparams}
\end{figure*}

\figref{fig:diffbloatmax} shows the reachable sets for the benchmark ISSF01
using different values of $\varepsilon_{\text{max}}{}$:
The smaller the defined error bound, the tighter are the reachable sets.
We also recognize the linear increase of the admissible error bound over time
as the computed sets differ more towards the end of the time horizon.
\figref{fig:ISSFparams} shows the evolution of all algorithm parameters{}
corresponding to the different values of $\varepsilon_{\text{max}}{}$:
A higher value for $\varepsilon_{\text{max}}{}$ yields a larger initial $\Delta t$
and a smaller total number of steps.
Note that the switching between previously-computed values of $\Delta t$ does not add
any computations since the sets are read from memory once they are computed.
The number of Taylor terms $\eta$ is chosen jointly with $\Delta t$
and increases towards the end of the time horizon to facilitate larger $\Delta t$.
We also observe that the zonotope order $\rho{}$
reaches a higher maximum for smaller values of $\varepsilon_{\text{max}}{}$
since in that case we cannot reduce as much as for larger $\varepsilon_{\text{max}}{}$.

\paragraph{Genetic Algorithm Comparison}
\label{par:GA}

Finally, we want to compare our approach to
a genetic algorithm searching for $\Delta t$, $\eta$, and~$\rho$.
To this end, we use the \emph{MATLAB} built-in genetic algorithm function.
While the parameters $\eta$ and $\rho$ are fixed,
we model the time step size by a polynomial up to order 2:
$\Delta t(t) = a + bt + ct^2$.
We restrict these parameters by the ranges
$\eta \in [1, 10]$, $\rho \in [2, 1000]$,
$a \in [0.0003,0.3], b \in [-0.1,0.1], c \in [-0.033,0.033]$.
The chosen bounds for $a,b$, and $c$ restrict $\Delta t$ to suitable curves.
Higher orders did not provide any benefits.

In order to establish a level playing field,
we terminate once the obtained reachable set is within the box enclosure
of the reachable set of the adaptive algorithm enlarged by 10\%.
For computational efficiency interval over-approximations were used for this comparison.
The cost function is chosen as by the maximum distance
to the enlarged adaptive reachable set over all dimensions.

The parameters specific to the genetic algorithm
have been set as follows:
We enable an infinite number of generations
with a maximum of 3 stall generations.
We aim to speed up the convergence by setting only 10 members per generation
as the evaluation of a single member is costly in higher dimensions.
For the members of the next generation, we use a standard crossover fraction of 0.75
and set the elite count to 1,
carrying the best solution over to the next generation.

We applied \algref{alg:full} and the genetic algorithm
on 50 randomly-generated systems per dimension.
\tabref{tab:GA_systems} compares the average computation time
over varying dimensions of \algref{alg:full}
to the time the genetic algorithms takes until convergence.
The results show that \algref{alg:full} outspeeds all genetic algorithms.
Since \algref{alg:full} tunes the algorithm parameters{} during runtime,
we only need a single iteration for the computation of the reachable set.
Contrary, the genetic algorithms run over many generations
repeatedly computing the reachable set while iteratively improving the solution
by means of the cost function.
This process is far more time-consuming than
the overhead caused by the adaptive parameter tuning.
The genetic algorithm using the constant polynomial for $\Delta t$
is faster than the higher-order polynomials as they re-compute
auxiliary reachable sets due to the non-constant time step size.

\vspace{-0.3cm}

\begin{table}[ht]
\begin{center}
	\caption{Computation time for \algref{alg:full} and the genetic algorithm (GA)
		averaged over 50 randomly-generated systems per dimension.}
	\label{tab:GA_systems}
	\begin{tabular}{l @{\hspace{20pt}}
	                c @{\hspace{8pt}}
	                c @{\hspace{8pt}}
	                c @{\hspace{8pt}}
	                c @{\hspace{8pt}}
	                c @{\hspace{8pt}}
	                c @{\hspace{8pt}}
	                c}
		\toprule
		& \multicolumn{7}{c}{\textbf{Dimension}} \\
		& 5 & 10 & 15 & 20 & 25 & 30 & 40 \\ \cmidrule{2-8}
		\algref{alg:full}            & 0.14s &  0.22s & 0.40s & 0.58s & 1.0s & 1.9s & 6.7s \\
		GA (order: 0) & 1.6s  &  3.4s  & 7.0s  & 10s   & 20s  & 30s  & 50s \\
		GA (order: 1) & 4.1s  &  9.4s  & 13s   & 21s   & 28s  & 40s  & 70s \\
		GA (order: 2) & 5.1s  &  10s   & 14s   & 21s   & 42s  & 58s  & 97s \\
		\bottomrule
	\end{tabular}
\end{center}
\vspace{-15pt}
\end{table}

\paragraph{Discussion}
\label{par:discussion}

Our framework can also be applied to
other computations of reachable sets and other set representations.
In order to guarantee convergence and termination for all $\varepsilon_{\text{max}}{}~\in~\mathbb{R}{}$,
\thmref{thm:alg:full} has to hold for the applied error terms,
similarly to \sectref{ssec:verass} for the presented implementation:
A tool developer has to modify
$\Phi_{\text{prop,max}}{}$ and $\Phi_{\text{set}}{}$ which are, e.g., the number of template directions
when using template polyhedra.
An corresponding error term for the set representation has to be defined.

The choice of $\errOp{\cdot}$ in \eqref{eq:bloat}
does not add much over-approximation as can be observed from
comparing the ranges for the time step size $[\Delta t_\text{min}, \Delta t_\text{max}]$
with the manually-tuned $\Delta t$ by \textit{CORA} in \tabref{tab:ARCH},
we see that \algref{alg:full} chooses a similar time step size,
yielding comparable results both in terms of the tightness
and the computational efficiency.

The only remaining parameter for the practitioner to set is the error $\varepsilon_{\text{max}}{}$.
As shown in \figref{fig:diffbloatmax}, increasing the value of $\varepsilon_{\text{max}}{}$
results in a more over-approximative reachable set and vice versa.
Thus, the setting of $\varepsilon_{\text{max}}{}$ is intuitive
and can easily be adjusted for any system.

\section{Preliminaries}
\label{sec:prelim}

To concisely explain the novelties of our paper,
we first introduce some preliminaries.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Notation}
\label{ssec:notation}

Vectors are denoted by lower-case letters
and matrices by upper-case letters.
Square matrices of zeros are denoted by
${\zeroMat{n} \in \mathbb{R}{n \times n}}$
and the identity matrix by ${\id{n} \in \mathbb{R}{n \times n}}$.
Given a vector ${a \in \mathbb{R}{n}}$, $a_i$ refers to the $i$-th entry.
An  $n$-dimensional interval is denoted by
$\mathcal{I}{} = [\ubar{l},\bar{l}] \subset \mathbb{R}{n}$,
where $\ubar{l}_i \leq \bar{l}_i$, $\forall i \in \{1,...,n\}$.
The operators $\inf(\mathcal{I}{}) = \ubar{l} \in \mathbb{R}{n}$
and $\sup(\mathcal{I}{}) = \bar{l} \in \mathbb{R}{n}$
return the infimum and the supremum of
an interval $\mathcal{I}{} = [\ubar{l},\bar{l}]$.
Interval matrices are denoted in boldface:
$\mathbf{I}=[\ubar{I},\bar{I}] \in \mathbb{R}{m \times n}$,
where $\inf(\mathbf{I}) = \ubar{I} \in \mathbb{R}{m \times n}$
and $\sup(\mathbf{I}) = \bar{I} \in \mathbb{R}{m \times n}$.
The linear map is written without any operator between the operands,
the Minkowski sum is denoted by $\oplus$,
and the convex hull operator by $\CH{\cdot}$.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Reachability Analysis of Linear Systems}
\label{ssec:ra_linsys}

The presented technique for adaptive self-parametrization
is applied to linear time-invariant (LTI) systems
with initial states bounded by $\mathcal{X}^0{} \subset \mathbb{R}{n}$
and inputs bounded by $\mathcal{U}{} \subset \mathbb{R}{n}$:
\begin{align}
	\begin{split}  \label{eq:LTI}
	&\dot{x}(t) = A x(t) + u(t) \; , \\
	&\text{with} \quad x(0) \in \mathcal{X}^0{} \subset \mathbb{R}{n}, \;
		\forall t: u(t) \in \mathcal{U}{} \subset \mathbb{R}{n} \; ,
	\end{split}
\end{align}
where $A \in \mathbb{R}{n \times n}$ is the system matrix,
$x(t) \in \mathbb{R}{n}$ is the state vector,
and $u(t) \in~\mathbb{R}{n}$ is the input vector.

The above system also encompasses systems
$\dot{x}(t) = A x(t) + B \tilde{u}(t)$,
where $B \in \mathbb{R}{m \times n}$ is the input matrix,
as we can set $\mathcal{U} = \{B\tilde{u} \, | \, \tilde{u} \in \mathcal{D} \subset \mathbb{R}{m}\}$.
Let us introduce $\xi \big(t; x_0, u(\cdot) \big)$
as the solution of \eqref{eq:LTI} to define the exact reachable set
$\Rexti{[0,T]}$ of \eqref{eq:LTI} over the time horizon $t \in [0,T]$:
\begin{align*}
	\Rexti{[0,T]} &= \Big\{ \xi \big(t; x_0, u(\cdot) \big) \, \Big| \, x_0 \in \mathcal{X}^0{}, \\
	&\qquad \quad \forall \tau \in [0,t]: u(\tau) \in \mathcal{U}{}, \, t \in [0,T] \Big\} \; .
\end{align*}

Due to the superposition principle, $\Rexti{[0,T]}$ is
the sum of the homogeneous solution~$\HRexti{[0,T]}$
and the inhomogeneous solution~$\PRexti{[0,T]}$ \cite[(3.6)]{Althoff2010a}:
\begin{equation*}
	\Rexti{[0,T]} = \HRexti{[0,T]} \oplus \PRexti{[0,T]} \; ,
\end{equation*}
where
\begin{align*}
	\HRexti{[0,T]} &= \Big\{ e^{A t} x_0 \, \Big|
		\, x_0 \in \mathcal{X}^0{}, \, t \in [0,T] \Big\} \; , \\
	\PRexti{T} &= \bigg\{ \int_0^T e^{A (T-\tau)} u(\tau) \mathrm{d}\tau \, \bigg|
			\, u(\tau) \in \mathcal{U}{} \bigg\} \; .
\end{align*}
If $0 \in \mathcal{U}{}$, $\PRexti{[0,T]} = \PRexti{T}$
and the extension for $0 \notin \mathcal{U}{}$
is shown in \cite[Sec. 3.2.2]{Althoff2010a}.
To limit the over-approximation,
the time horizon $[0,T]$ is discretized by $K$ time steps
$\Delta t_i = t_{i+1} - t_i > 0$, $i \in \{0,...,K-1\}$,
where $T = \sum_{i=0}^{K-1} \Delta t_i$
so that
$\bigcup_{i=0}^{K-1}~\Rti{[t_i,t_{i+1}]}$.
We compute the reachable sets from the first time interval
\cite[Chap. 4]{LeG2009}:
\begin{align}
	\begin{split} \label{eq:HRPR_start}
		\HRti{[t_i,t_{i+1}]} &= e^{A t_i} \HRti{[0,\Delta t_i]} \; ,  \\
		\PRti{[0,t_{i+1}]} &= \PRti{[0, t_i]} \oplus
			\underbrace{e^{A t_i} \PRti{[0,\Delta t_i]}}_%
			{=: \, \PRtiunderbrace{[t_i,t_{i+1}]}} \; .
	\end{split}
\end{align}
This method adapts seamlessly to varying time step sizes,
making it the preferred choice for adaptive parameter tuning
compared to methods propagating the reachable set from the previous time step.

\section{Self-Parametrization}
\label{sec:selfparam}

In this section, we introduce our novel algorithm for adaptive parameter tuning.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Overview of Algorithm Parameters{}}
\label{ssec:paramoverview}

The reachability algorithm from \sectref{ssec:ra_linsys}
depends on different algorithm parameters{},
which we divide into:
\begin{itemize}
	\item \textbf{Time Step Size:}
		The propagation in~\eqref{eq:HRPR_start} can be applied to any series
		of $\Delta t_i$ adding up to $T$.
		Large $\Delta t$ speed up the computation, while small $\Delta t$ increase the precision.
	\item \textbf{Propagation Parameters{}:}
		The tightness of the computed reachable sets
		also depends on parameters approximating the dynamics $\Phi_{\text{prop}}{}$,
		e.g., determining the precision of computing $e^{At}$.
		We introduce the operator $\incrOp{\Phi_{\text{prop}}{}}$
		which increments the parameters $\Phi_{\text{prop}}{}$ towards values
		which result in a tighter reachable set,
		as well as the operator $\resetOp{\Phi_{\text{prop}}{}}$
		which resets~$\Phi_{\text{prop}}{}$ to the coarsest setting.
	\item \textbf{Set Representation:}
		We denote the number of scalar values needed to describe a set by~$\Phi_{\text{set}}{}$.
		Let us introduce the operator $\decrOp{\Phi_{\text{set}}{}}$ which decreases $\Phi_{\text{set}}{}$.
		If $\Phi_{\text{set}}{}$ has reached its minimum, it returns $\decrOp{\Phi_{\text{set}}{}} = \Phi_{\text{set}}{}$.
		We denote the over-approximation of a set $\mathcal{S}{}$
		by reducing its number of parameters
		to a desired value $\Phi_{\text{set}}{}$ by the operator $\redOp{\mathcal{S}{}}{\Phi_{\text{set}}{}}$.
\end{itemize}

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Error measures}
\label{ssec:errors}

As we can only compute over-approximations,
we are interested in measuring the over-approximation error.
Let the reachable set $\mathcal{H}^\mathcal{R}{}$ consist of summands, e.g.,
$\mathcal{H}^\mathcal{R}{} = \mathcal{H}^\mathcal{R}{}^{(1)} \oplus ... \oplus \mathcal{H}^\mathcal{R}{}^{(\nu)}$,
as shown for a concrete implementation in \sectref{ssec:setprop}.
This assumption is valid for almost all approaches, see e.g.
\cite{Althoff2010a,Girard2005,Girard2008b,FreLeG2011,BFFPSV2018}.
We sum over all indices $\mathcal{E}$ representing exactly-computed terms
to obtain $\mathcal{H}^\mathcal{R}_={} = \bigoplus_{i \in \mathcal{E}}^{\nu} \mathcal{H}^\mathcal{R}{}^{(i)}$
and the remaining terms are summed to $\mathcal{H}^\mathcal{R}_+{} = \bigoplus_{i \notin \mathcal{E}}^{\nu} \mathcal{H}^\mathcal{R}{}^{(i)}$.
The same is done for $\mathcal{P}^\mathcal{R}{}$, so that we obtain
\begin{align}
	\HRti{[0,\Delta t]} &= \HReqti{[0,\Delta t]} \oplus \HRplusti{[0,\Delta t]} \; ,
		\label{eq:HRtieqplus} \\
	\PRti{[0,\Delta t]} &= \PReqti{[0,\Delta t]} \oplus \PRplusti{[0,\Delta t]} \; .
		\label{eq:PRtieqplus}
\end{align}
We use this separation to realize a computationally efficient over-approximation
of the Hausdorff distance
\begin{align*}
	d_H (\mathcal{S}{}_1, \mathcal{S}{}_2)
		= \max &\Big\{ \sup_{x_1 \in \mathcal{S}{}_1}
				\big( \, \inf_{x_2 \in \mathcal{S}{}_2} \, \norm{x_1 - x_2} \,
			\big), \\
		&\quad \sup_{x_2 \in \mathcal{S}{}_2}
			\big( \, \inf_{x_1 \in \mathcal{S}{}_1} \norm{x_1 - x_2} \,
			\big) \Big\} \; .
\end{align*}
\begin{proposition} \label{prop:bloat}
Let ${\mathcal{S}_{=}}{} \subset \mathbb{R}{n}$ and $\mathcal{S}_{+}{} \subset \mathbb{R}{n}$
with ${0 \in \mathcal{S}_{+}{}}$ be non-empty compact, convex sets.
The Hausdorff distance between ${\mathcal{S}_{=}}{}$ and ${\mathcal{S}_{\text{tot}}}{} := {\mathcal{S}_{=}}{} \oplus \mathcal{S}_{+}{}$
can be over-approximated by
\begin{equation} \label{eq:bloat}
	d_H({\mathcal{S}_{=}}{},{\mathcal{S}_{\text{tot}}}{}) \leq \errOp{\mathcal{S}_{+}{}}
		:= r \big( \boxOp{\mathcal{S}_{+}{}} \big) \; ,
\end{equation}
where $\boxOp{\mathcal{S}_{+}{}}$ is the box enclosure of $\mathcal{S}_{+}{}$
and $r(\cdot)$ returns the radius of the smallest hypersphere
centered at the origin enclosing its argument.
\end{proposition}
\begin{proof}
Since ${\mathcal{S}_{\text{tot}}}{}$ encloses ${\mathcal{S}_{=}}{}$, we have
\begin{equation} \label{eq:dH_r}
	d_H({\mathcal{S}_{=}}{},{\mathcal{S}_{\text{tot}}}{})
		= \sup_{y \in {\mathcal{S}_{\text{tot}}}{}} \big( \, \inf_{x \in {\mathcal{S}_{=}}{}} \, \norm{x - y} \big) \; .
\end{equation}
We also know that the difference between ${\mathcal{S}_{\text{tot}}}{}$ and ${\mathcal{S}_{=}}{}$ is given by $\mathcal{S}_{+}{}$,
so that \eqref{eq:dH_r} can be equivalently written as
\begin{equation*}
	\sup_{y \in {\mathcal{S}_{\text{tot}}}{}} \big( \, \inf_{x \in {\mathcal{S}_{=}}{}} \norm{x - y} \, \big)
		\overset{0 \in \mathcal{S}_{+}{}}{=} \sup_{x \in \mathcal{S}_{+}{}} \norm{x} = r ( \mathcal{S}_{+}{} ) \; .
\end{equation*}
Obviously, $r(\mathcal{S}_{+}{}) \leq r\big(\boxOp{\mathcal{S}_{+}{}}\big)$
and therefore $d_H({\mathcal{S}_{=}}{},{\mathcal{S}_{\text{tot}}}{}) \leq \errOp{\mathcal{S}_{+}{}}$.
\end{proof}

Let us introduce user-defined upper bounds for errors:
$\varepsilon_{\mathcal{H},\text{max}}{}$ for the error of the homogeneous solution $\mathcal{H}^\mathcal{R}{}$ and
$\varepsilon_{\mathcal{P},\text{max}}{}$ for the error of the inhomogeneous solution $\mathcal{P}^\mathcal{R}{}$.
Our approach ensures that these errors are not exceeded.

We now compute the errors of the $i$-th time step.
For the error $\bloatHRti{[t_i,t_{i+1}]}$ and $\bloatPRti{[t_i,t_{i+1}]}$
we use the split in~\eqref{eq:HRtieqplus} and apply the error measure from \propref{prop:bloat}:
\begin{align}
	\bloatHRti{[t_i,t_{i+1}]} &:= \errOp{\HRplusti{[t_i,t_{i+1}]}} \label{eq:bloatHRti} \\
		&\geq d_H \big( \HReqti{[t_i,t_{i+1}]},\HRti{[t_i,t_{i+1}]} \big) \; , \nonumber \\
	\bloatPRti{[t_i,t_{i+1}]} &:= \errOp{\PRplusti{[t_i,t_{i+1}]}} \label{eq:bloatPRti} \\
		&\geq d_H \big( \PReqti{[t_i,t_{i+1}]},\PRti{[t_i,t_{i+1}]} \big) \; . \nonumber
\end{align}
From \eqref{eq:HRPR_start}, we see that the computation of $\varepsilon_{\mathcal{H}}{}$ does not depend on previous time intervals.
Contrary, the inhomogeneous solution $\PRti{[0,t_i]}$
accumulates over time, see \eqref{eq:HRPR_start},
and with it the error $\varepsilon_{\mathcal{P}}{}$:
\begin{equation}
	\bloatPRti{[0,t_{i+1}]} = \bloatPRti{[0,t_i]} + \bloatPRti{[t_i,t_{i+1}]} \; ,
		\label{eq:bloatPRti_acc}
\end{equation}
where $\bloatPRti{[0,0]} = 0$.
Let us introduce a mild assumption for the errors $\varepsilon_{\mathcal{H}}{}$ and $\varepsilon_{\mathcal{P}}{}$,
which will be justified in \sectref{ssec:verass}.

\begin{assumption} \label{ass:bloatHRPR} (Convergence of $\varepsilon_{\mathcal{H}}{}$ and $\varepsilon_{\mathcal{P}}{}$)
We neglect floating-point errors so that $\bloatHRti{[t_i,t_{i+1}]} \rightarrow 0$
and $\bloatPRti{[t_i,t_{i+1}]} \rightarrow 0$ for $\Delta t_i \rightarrow 0$. \hfill $\square$
\end{assumption}
\assref{ass:bloatHRPR} is mild, since it practically holds
for all other current approaches, see e.g. \cite{Althoff2015a, FreLeG2011, BogFor2019}.
For further derivations we need the following definition.
\begin{definition} \label{def:superlinear} (Superlinear decrease)
	The set-based evaluation $f(\mathcal{S}{}(t)) = \{ f(x) \, | \, x \in \mathcal{S}{}(t) \}$
	is a continuous function $f: \mathbb{R}{n} \rightarrow \mathbb{R}{}$
	decreases superlinearly if
	\begin{equation*}
		\hspace{8pt} \forall \varphi \in (0,1):
			f\big(\mathcal{S}{}\big([t,t + \varphi \Delta t]\big)\big)
				\subseteq \varphi \, f\big(\mathcal{S}{}\big([t,t + \Delta t]\big)\big) \; .
		\hspace{5pt} \square
	\end{equation*}
\end{definition}
The following proposition addresses satisfying $\varepsilon_{\mathcal{P},\text{max}}{}$.

\begin{proposition} \label{prop:bloatPRtiadm}
Let the error $\bloatPRti{[t_i,t_{i+1}]}$
decrease superlinearly according to \defref{def:superlinear}.
For any given $\varepsilon_{\mathcal{P},\text{max}}{} > 0$, there exists a sequence of
time intervals $[t_i,t_{i+1}]$ so that
\begin{equation*}
	\sum_i \bloatPRti{[t_i,t_{i+1}]} \leq \varepsilon_{\mathcal{P},\text{max}}{} \; .
\end{equation*}
\end{proposition}

\noindent \textit{Proof.}
We define an admissible error $\bloatPRtiadm{i}$ for each step:
\begin{equation} \label{eq:bloatPRtiadm}
	\bloatPRti{[t_i,t_{i+1}]} \leq
	\bloatPRtiadm{i} := \frac{\varepsilon_{\mathcal{P},\text{max}}{} - \bloatPRti{[0,t_i]}}{T - t_i} \, \Delta t_i .
\end{equation}
For $\Delta t_i \rightarrow 0$, $\bloatPRtiadm{i}$ converges to 0,
as does $\bloatPRti{[t_i,t_{i+1}]}$ by \assref{ass:bloatHRPR}.
Moreover, $\bloatPRtiadm{i}$ decreases linearly in~$\Delta t_i$,
whereas $\bloatPRti{[t_i,t_{i+1}]}$ decreases superlinearly by assumption.
Hence, we can always obtain a $\Delta t_i$ so that the inequality in \eqref{eq:bloatPRtiadm} holds.
We sum each side in \eqref{eq:bloatPRtiadm} to
\begin{equation} \label{eq:bloatPRsumti}
	\bloatPRti{[0,t_i]} \leq \sum_{j=0}^{i-1} \bloatPRtiadm{j}
\end{equation}
which is subsequently used to bound $\bloatPRtiadm{i}$ and $\varepsilon_{\mathcal{P},\text{max}}{}$:
\begin{align*}
	\hspace{26pt}
	\bloatPRtiadm{i}
		&\overset{\eqref{eq:bloatPRtiadm}}{=} \big( \varepsilon_{\mathcal{P},\text{max}}{} - \bloatPRti{[0,t_i]} \big) \frac{\Delta t_i}{T - t_i} \\
		&\overset{\eqref{eq:bloatPRsumti}}{\leq} \big( \varepsilon_{\mathcal{P},\text{max}}{} - \sum_{j=0}^{i-1} \bloatPRtiadm{j} \big) \frac{\Delta t_i}{T - t_i} \\
	\Rightarrow \varepsilon_{\mathcal{P},\text{max}}{}
		&\geq \sum_{j=0}^{i-1} \bloatPRtiadm{j} + \underbrace{\frac{T - t_i}{\Delta t_i}}_{\geq 1} \bloatPRtiadm{i} \\
		&\geq \sum_{j=0}^{i} \bloatPRtiadm{j} \geq \sum_{i} \bloatPRti{[t_i, t_{i+1}]} \; .
	\hspace{26pt}
	\square
\end{align*}

The Minkowski sum in~\eqref{eq:HRPR_start} typically increases the representation size of the resulting set.
To counteract this growth, we enclose~$\mathcal{P}^\mathcal{R}{}$ in~\eqref{eq:HRPR_start}
by a set specified by less parameters which over-approximates the original set
by the error $\varepsilon_{\mathcal{S}}{}$ accumulating
as the error $\varepsilon_{\mathcal{P}}{}$ of $\mathcal{P}^\mathcal{R}{}$:
\begin{align}
	&\bloatSetti{[0,t_{i+1}]} = \bloatSetti{[0,t_i]} + \bloatSetti{[t_i,t_{i+1}]} \; ,
		\label{eq:bloatSetti_prop} \\
	&\text{with} \quad \bloatSetti{[0,0]} = 0 \; . \nonumber
\end{align}
Similarly to \propref{prop:bloatPRtiadm},
we define an admissible error $\bloatSettiadm{i}$:
\begin{equation} \label{eq:bloatSettiadm}
	\bloatSetti{[t_i,t_{i+1}]} \leq \bloatSettiadm{i} := \frac{\varepsilon_{\mathcal{S},\text{max}}{} - \bloatSetti{[0,t_i]}}{T - t_i} \, \Delta t_i \; .
\end{equation}
If we do not over-approximate the set representation
in the $i$-th time step, we have $\bloatSetti{[t_i,t_{i+1}]} = 0$.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{Adaptive Parameter Tuning}
\label{ssec:tuning}

We now present our generic framework for adaptive parameter tuning
in \algref{alg:full}.
The two subroutines, \algref{alg:postTSTT} and \algref{alg:propInhom},
will be presented thereafter.
First, we initialize $\Delta t$ and the scaling factor $\mu$ (\lineref{alg:full:inits}),
as well as the accumulating errors $\varepsilon_{\mathcal{P}}{}$ and $\varepsilon_{\mathcal{S}}{}$
(\lineref{alg:full:initbloat})
for the first iteration of the main loop
(\linerangeref{alg:full:mainloop_start}{alg:full:mainloop_end}).
In every iteration, we first obtain the parameters
$\Delta t_i$ and $\AppModPidx{i}$ (\lineref{alg:full:postTSTT})
from \algref{alg:postTSTT}.
Then, we calculate the reachable sets
$\HRti{[0,\Delta t_i]}$ and $\PRti{[0,\Delta t_i]}$ (\lineref{alg:full:HRPR0Deltati}),
which are subsequently used to obtain $\HRti{[t_i,t_{i+1}]}$
and $\PRti{[0,t_{i+1}]}$ (\lineref{alg:full:HRPR_prop}).
In \algref{alg:propInhom}, we obtain the parameters $\SetReprPidx{i}$,
which are used to recalculate $\PRti{[0,t_{i+1}]}$ (\lineref{alg:full:propInhom}).
At the end of each step, we add the homogeneous and inhomogeneous solution
to obtain the reachable set $\Rti{[t_i, t_{i+1}]}$ (\lineref{alg:full:Rti}).
Finally, after the calculation of the reachable set for each time interval,
the reachable set of the whole time horizon $\Rti{[0,T]}$
is obtained by unifying partial sets (\lineref{alg:full:R0T}).

\begin{theorem} \label{thm:alg:full}
\algref{alg:full} terminates and the overall error
stays below a user-defined error bound $\varepsilon_{\text{max}}{} \in \mathbb{R}{+}$,
assuming the reduction of the set representation in each step is optional.
\end{theorem}

\FloatBarrier

\begin{algorithm}[h!tb]
	\caption{Fully automated parameter tuning} \label{alg:full}
	\textbf{Input:} $\varepsilon_{\mathcal{H},\text{max}}{}, \varepsilon_{\mathcal{P},\text{max}}{}, \varepsilon_{\mathcal{S},\text{max}}{}, T$

	\textbf{Output:} $\Rti{[0,T]}$
	\begin{algorithmic}[1]
		\State $t \gets 0, i \gets 0$
		\State $\mu{} \gets 0.9, \Delta t_{-1} \gets T \, \mu{}$
			\label{alg:full:inits}
		\State $\bloatPRti{[0,0]} \gets 0$, $\bloatSetti{[0,0]} \gets 0$
			\label{alg:full:initbloat}
		\While{$t < T$}
			\label{alg:full:mainloop_start}
			\State $\Delta t_i, \AppModPidx{i}, \bloatPRti{[0,t_{i+1}]} \gets$ \algref{alg:postTSTT}
				\label{alg:full:postTSTT}
			\State calc. $\HRti{[0,\Delta t_i]}, \PRti{[0,\Delta t_i]}$ using $\Delta t_i$, $\AppModPidx{i}$
				\label{alg:full:HRPR0Deltati}
			\State calc. $\HRti{[t_i,t_{i+1}]}, \PRti{[0,t_{i+1}]}$ acc. to \eqref{eq:HRPR_start}
				\label{alg:full:HRPR_prop}
			\State $\PRti{[0,t_{i+1}]}, \bloatSetti{[0,t_{i+1}]} \gets$ \algref{alg:propInhom}
				\label{alg:full:propInhom}
			\State $\Rti{[t_i, t_{i+1}]} = \HRti{[t_i,t_{i+1}]} \oplus \PRti{[0,t_{i+1}]}$
				\label{alg:full:Rti}
			\State $t \gets t + \Delta t_i$, $i \gets i+1$
				\label{alg:full:increaset}
		\EndWhile \label{alg:full:mainloop_end}
		\State $\Rti{[0,T]} \gets \bigcup_{j=0}^{i-1} \Rti{[t_j, t_{j+1}]}$
			\label{alg:full:R0T}
	\end{algorithmic}
\end{algorithm}

\begin{proof}
We first show that the algorithm terminates.
The while-loop in \algref{alg:postTSTT} always terminates if
\begin{equation*}
	\forall i: \bloatPRti{[t_i, t_{i+1}]} \leq \bloatPRtiadm{i}
	\; , \; \bloatHRti{[t_i, t_{i+1}]} \leq \varepsilon_{\mathcal{H},\text{max}}{}	\; .
\end{equation*}
The first inequality can be satisfied as shown in the proof of \propref{prop:bloatPRtiadm}.
From \lineref{alg:postTSTT:decreaseDeltat} in \algref{alg:postTSTT}
and under \assref{ass:bloatHRPR}, we have
$\forall i: \Delta t_i \rightarrow 0 \Rightarrow \bloatHRti{[t_i,t_{i+1}]} \rightarrow 0$.
Thus, every bound $\varepsilon_{\mathcal{H},\text{max}}{}$ is satisfiable.
The while-loop in \algref{alg:propInhom} terminates as we will
either exceed the admissible error bound $\bloatSettiadm{i}$
during the continuous reduction or exit the loop
if the set parameters $\Phi_{\text{set}}{}$ have been reduced as much as possible.
Finally, the main loop (\linerangeref{alg:full:mainloop_start}{alg:full:mainloop_end})
terminates as the time~$t$ monotonically increases (\lineref{alg:full:increaset}),
eventually reaching~$T$.

To prove that the error stays below a user-defined bound,
we split the error $\varepsilon_{\text{max}}{}$
in $\varepsilon_{\text{max}}{} = \varepsilon_{\mathcal{H},\text{max}}{} + \varepsilon_{\mathcal{P},\text{max}}{} + \varepsilon_{\mathcal{S},\text{max}}{}$.
All individual bounds
are satisfiable:
We have already shown the satisfiability of any $\varepsilon_{\mathcal{H},\text{max}}{}$
in the beginning of this proof.
\propref{prop:bloatPRtiadm} shows that any limit $\varepsilon_{\mathcal{P},\text{max}}{}$ is satisfiable.
Lastly, any bound~$\varepsilon_{\mathcal{S},\text{max}}{}$ can be satisfied since
we can choose to not over-approximate the set representation.
\end{proof}

\noindent We now present two algorithms performing the adaptive parameter tuning:
As an overview, \algref{alg:postTSTT} adapts
the parameters $\Delta t_i$ and $\AppModPidx{i}$,
while \algref{alg:propInhom} adapts the parameters $\SetReprPidx{i}$
concurrently to the propagation of $\mathcal{P}^\mathcal{R}{}$.

We first present \algref{alg:postTSTT}:
As the main idea, every $\Phi_{\text{prop}}{}$ is checked for an iteratively decreasing $\Delta t_i$
until the error bounds $\varepsilon_{\mathcal{H},\text{max}}{}$ and $\bloatPRtiadm{i}$ are satisfied.
Initially, we increase the time step size by division with the factor $\mu{}$
(\lineref{alg:postTSTT:startDeltat}).
Furthermore, we reset~$\AppModPidx{i}$ to its coarsest setting (\lineref{alg:postTSTT:startDeltat})
and calculate the admissible error $\bloatPRtiadm{i}$ (\lineref{alg:postTSTT:bloatPRtiadm}).
In the loop (\linerangeref{alg:postTSTT:loop_start}{alg:postTSTT:loop_end}),
we increment $\AppModPidx{i}$ (\lineref{alg:postTSTT:incrAppModP})
and calculate the errors $\bloatHRti{[t_i,t_{i+1}]}$ and $\bloatPRti{[t_i,t_{i+1}]}$
based on the error sets $\HRplusti{[t_i,t_{i+1}]}$ and $\PRplusti{[t_i,t_{i+1}]}$
(\linerangeref{alg:postTSTT:errors_start}{alg:postTSTT:errors_end}).
These errors are then compared to their respective error bounds
(\lineref{alg:postTSTT:loop_end}).
If $\AppModPidx{i}$ reaches $\Phi_{\text{prop,max}}{}$ (\lineref{alg:postTSTT:AppModPmax}),
we decrease $\Delta t_i$, reset $\AppModPidx{i}$ (\lineref{alg:postTSTT:decreaseDeltat}),
and recalculate the admissible bound $\bloatPRtiadm{i}$
(\lineref{alg:postTSTT:updatebloatPRtiadm})
before restarting the loop (\lineref{alg:postTSTT:continue}).
After the loop is finished, the accumulated error $\bloatPRti{[0,t_{i+1}]}$
is computed (\lineref{alg:postTSTT:bloatPRti_acc}).

\FloatBarrier

\begin{algorithm}[h!tb]
	\caption{Adapt values for $\Delta t, \Phi_{\text{prop}}{}$} \label{alg:postTSTT}
	\textbf{Input:} $\varepsilon_{\mathcal{H},\text{max}}{}, \varepsilon_{\mathcal{P},\text{max}}{}, \bloatPRti{[0,t_i]}, \Delta t_{i-1}, \mu{}, T$

	\textbf{Output:} $\Delta t_i, \AppModPidx{i}, \bloatPRti{[0,t_{i+1}]}$
	\begin{algorithmic}[1]
		\State $\Delta t_i \gets \frac{\Delta t_{i-1}}{\mu{}}$, $\resetOp{\AppModPidx{i}}$
			\label{alg:postTSTT:startDeltat}
		\State calc. $\bloatPRtiadm{i}$ acc. to \eqref{eq:bloatPRtiadm}
			\label{alg:postTSTT:bloatPRtiadm}
		\State \textbf{do}
			\label{alg:postTSTT:loop_start}
			\State $\quad$ $\incrOp{\AppModPidx{i}}$
				\label{alg:postTSTT:incrAppModP}
			\State $\quad$ \textbf{if} $\Phi_{\text{prop}}{} \geq \Phi_{\text{prop,max}}{}$
				\label{alg:postTSTT:AppModPmax}
				\State $\qquad$ $\Delta t_i \gets \Delta t_i \mu{}$, $\resetOp{\AppModPidx{i}}$
					\label{alg:postTSTT:decreaseDeltat}
				\State $\qquad$ calc. $\bloatPRtiadm{i}$ acc. to \eqref{eq:bloatPRtiadm}
					\label{alg:postTSTT:updatebloatPRtiadm}
				\State $\qquad$ \textbf{continue}
					\label{alg:postTSTT:continue}
			\State $\quad$ \textbf{end if}
			\State $\quad$ calc. $\HRplusti{[t_i,t_{i+1}]}, \bloatHRti{[t_i, t_{i+1}]}$
				\label{alg:postTSTT:errors_start}
			\State $\quad$ calc. $\PRplusti{[t_i,t_{i+1}]}, \bloatPRti{[t_i, t_{i+1}]}$
				\label{alg:postTSTT:errors_end}
		\State \textbf{while} $\bloatHRti{[t_i, t_{i+1}]} > \varepsilon_{\mathcal{H},\text{max}}{}
			\wedge \bloatPRti{[t_i, t_{i+1}]} > \bloatPRtiadm{i}$
				\label{alg:postTSTT:loop_end}
		\State calc. $\bloatPRti{[0,t_{i+1}]}$ acc. to \eqref{eq:bloatPRti_acc}
			\label{alg:postTSTT:bloatPRti_acc}
	\end{algorithmic}
\end{algorithm}

\FloatBarrier

The adaptation of~$\Phi_{\text{set}}{}$ is shown in \algref{alg:propInhom}:
In order to decrease the computation time, the number of parameters
describing $\mathcal{P}^\mathcal{R}{}$ is iteratively reduced from a high-precision representation
towards lower precision.
First, we calculate the admissible error~$\bloatSettiadm{i}$
(\lineref{alg:propInhom:bloatSettiadm}) for the given~$\Delta t_i$.
The operator $\paramOp{\Phi_{\text{set}}{}}$ extracts
the number of stored parameters~$\SetReprPidx{i}$
from $\PRti{[0,t_{i+1}]}$ (\lineref{alg:propInhom:readSetReprP}).
Inside the loop (\linerangeref{alg:propInhom:loop_start}{alg:propInhom:loop_end}),
we iteratively decrease $\SetReprPidx{i}$ (\lineref{alg:propInhom:decrSetReprP})
and over-approximate $\PRti{[0,t_{i+1}]}$
by reducing the number of stored parameters to $\SetReprPidx{i}$
(\lineref{alg:propInhom:reduce}).
From the difference between the original set and the set over-approximated by reduction,
we calculate the induced error~$\bloatSetti{[t_i,t_{i+1}]}$ (\lineref{alg:propInhom:bloatSetti})
and compare it to the admissible error (\lineref{alg:propInhom:loop_end}).
After the loop,
the accumulated error for the set representation $\bloatSetti{[0,t_{i+1}]}$
is computed (\lineref{alg:propInhom:bloatSetti_prop}).

\FloatBarrier

\begin{algorithm}[h!tb]
	\caption{Adapt values for $\Phi_{\text{set}}{}$} \label{alg:propInhom}
	\textbf{Input:} $\varepsilon_{\mathcal{S},\text{max}}{}, \bloatSetti{[0,t_i]}, \PRti{[0,t_{i+1}]}, \Delta t_i, T$

	\textbf{Output:} $\PRti{[0,t_{i+1}]}, \bloatSetti{[0,t_{i+1}]}$
	\begin{algorithmic}[1]
		\State calc. $\bloatSettiadm{i}$ acc. to \eqref{eq:bloatSettiadm}
			\label{alg:propInhom:bloatSettiadm}
		\State $\SetReprPidx{i} \gets \paramOp{\PRti{[0,t_{i+1}]}}$
			\label{alg:propInhom:readSetReprP}
		\State \textbf{do}
			\label{alg:propInhom:loop_start}
			\State $\quad$ $\decrOp{\SetReprPidx{i}}$
				\label{alg:propInhom:decrSetReprP}
			\State $\quad$ $\redOp{\PRti{[0,t_{i+1}]}}{\SetReprPidx{i}}$
				\label{alg:propInhom:reduce}
			\State $\quad$ calc. $\bloatSetti{[t_i,t_{i+1}]}$
				\label{alg:propInhom:bloatSetti}
		\State \textbf{while} $\bloatSetti{[t_i,t_{i+1}]} < \bloatSettiadm{i}
			\vee \decrOp{\Phi_{\text{set}}{}} = \Phi_{\text{set}}{}$
			\label{alg:propInhom:loop_end}
		\State calc. $\bloatSetti{[0,t_{i+1}]}$ acc. to \eqref{eq:bloatSetti_prop}
			\label{alg:propInhom:bloatSetti_prop}
	\end{algorithmic}
\end{algorithm}

\FloatBarrier

The presented framework for adaptive parameter tuning
can be applied to any reachability algorithm.
In the next section, we will provide an example implementation
and prove the assumptions underlying \thmref{thm:alg:full}.

\@startsection{subsection}{2}{\z@}{-3.25ex\@plus -1ex \@minus -.2ex}{1.5ex \@plus .2ex{%

redpajama_set_name: RedPajamaArXiv