**Dana Fisman Grigore Rosu (Eds.)**

# **Tools and Algorithms for the Construction and Analysis of Systems**

**28th International Conference, TACAS 2022 Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2022 Munich, Germany, April 2–7, 2022 Proceedings, Part II**

# Lecture Notes in Computer Science 13244

Founding Editors

Gerhard Goos, Germany Juris Hartmanis, USA

# Editorial Board Members

Elisa Bertino, USA Wen Gao, China Bernhard Steffen , Germany Gerhard Woeginger , Germany Moti Yung , USA

# Advanced Research in Computing and Software Science Subline of Lecture Notes in Computer Science

Subline Series Editors

Giorgio Ausiello, University of Rome 'La Sapienza', Italy Vladimiro Sassone, University of Southampton, UK

Subline Advisory Board

Susanne Albers, TU Munich, Germany Benjamin C. Pierce, University of Pennsylvania, USA Bernhard Steffen , University of Dortmund, Germany Deng Xiaotie, Peking University, Beijing, China Jeannette M. Wing, Microsoft Research, Redmond, WA, USA More information about this series at https://link.springer.com/bookseries/558

# Tools and Algorithms for the Construction and Analysis of Systems

28th International Conference, TACAS 2022 Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2022 Munich, Germany, April 2–7, 2022 Proceedings, Part II

Editors Dana Fisman Ben-Gurion University of the Negev Be'er Sheva, Israel

Grigore Rosu University of Illinois Urbana-Champaign Urbana, IL, USA

ISSN 0302-9743 ISSN 1611-3349 (electronic) Lecture Notes in Computer Science ISBN 978-3-030-99526-3 ISBN 978-3-030-99527-0 (eBook) https://doi.org/10.1007/978-3-030-99527-0

© The Editor(s) (if applicable) and The Author(s) 2022. This book is an open access publication.

Open Access This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

# ETAPS Foreword

Welcome to the 25th ETAPS! ETAPS 2022 took place in Munich, the beautiful capital of Bavaria, in Germany.

ETAPS 2022 is the 25th instance of the European Joint Conferences on Theory and Practice of Software. ETAPS is an annual federated conference established in 1998, and consists of four conferences: ESOP, FASE, FoSSaCS, and TACAS. Each conference has its own Program Committee (PC) and its own Steering Committee (SC). The conferences cover various aspects of software systems, ranging from theoretical computer science to foundations of programming languages, analysis tools, and formal approaches to software engineering. Organizing these conferences in a coherent, highly synchronized conference program enables researchers to participate in an exciting event, having the possibility to meet many colleagues working in different directions in the field, and to easily attend talks of different conferences. On the weekend before the main conference, numerous satellite workshops took place that attract many researchers from all over the globe.

ETAPS 2022 received 362 submissions in total, 111 of which were accepted, yielding an overall acceptance rate of 30.7%. I thank all the authors for their interest in ETAPS, all the reviewers for their reviewing efforts, the PC members for their contributions, and in particular the PC (co-)chairs for their hard work in running this entire intensive process. Last but not least, my congratulations to all authors of the accepted papers!

ETAPS 2022 featured the unifying invited speakers Alexandra Silva (University College London, UK, and Cornell University, USA) and Tomáš Vojnar (Brno University of Technology, Czech Republic) and the conference-specific invited speakers Nathalie Bertrand (Inria Rennes, France) for FoSSaCS and Lenore Zuck (University of Illinois at Chicago, USA) for TACAS. Invited tutorials were provided by Stacey Jeffery (CWI and QuSoft, The Netherlands) on quantum computing and Nicholas Lane (University of Cambridge and Samsung AI Lab, UK) on federated learning.

As this event was the 25th edition of ETAPS, part of the program was a special celebration where we looked back on the achievements of ETAPS and its constituting conferences in the past, but we also looked into the future, and discussed the challenges ahead for research in software science. This edition also reinstated the ETAPS mentoring workshop for PhD students.

ETAPS 2022 took place in Munich, Germany, and was organized jointly by the Technical University of Munich (TUM) and the LMU Munich. The former was founded in 1868, and the latter in 1472 as the 6th oldest German university still running today. Together, they have 100,000 enrolled students, regularly rank among the top 100 universities worldwide (with TUM's computer-science department ranked #1 in the European Union), and their researchers and alumni include 60 Nobel laureates. The local organization team consisted of Jan Křetínský (general chair), Dirk Beyer (general, financial, and workshop chair), Julia Eisentraut (organization chair), and Alexandros Evangelidis (local proceedings chair).

ETAPS 2022 was further supported by the following associations and societies: ETAPS e.V., EATCS (European Association for Theoretical Computer Science), EAPLS (European Association for Programming Languages and Systems), and EASST (European Association of Software Science and Technology).

The ETAPS Steering Committee consists of an Executive Board, and representatives of the individual ETAPS conferences, as well as representatives of EATCS, EAPLS, and EASST. The Executive Board consists of Holger Hermanns (Saarbrücken), Marieke Huisman (Twente, chair), Jan Kofroň (Prague), Barbara König (Duisburg), Thomas Noll (Aachen), Caterina Urban (Paris), Tarmo Uustalu (Reykjavik and Tallinn), and Lenore Zuck (Chicago).

Other members of the Steering Committee are Patricia Bouyer (Paris), Einar Broch Johnsen (Oslo), Dana Fisman (Be'er Sheva), Reiko Heckel (Leicester), Joost-Pieter Katoen (Aachen and Twente), Fabrice Kordon (Paris), Jan Křetínský (Munich), Orna Kupferman (Jerusalem), Leen Lambers (Cottbus), Tiziana Margaria (Limerick), Andrew M. Pitts (Cambridge), Elizabeth Polgreen (Edinburgh), Grigore Roşu (Illinois), Peter Ryan (Luxembourg), Sriram Sankaranarayanan (Boulder), Don Sannella (Edinburgh), Lutz Schröder (Erlangen), Ilya Sergey (Singapore), Natasha Sharygina (Lugano), Pawel Sobocinski (Tallinn), Peter Thiemann (Freiburg), Sebastián Uchitel (London and Buenos Aires), Jan Vitek (Prague), Andrzej Wasowski (Copenhagen), Thomas Wies (New York), Anton Wijs (Eindhoven), and Manuel Wimmer (Linz).

I'd like to take this opportunity to thank all authors, attendees, organizers of the satellite workshops, and Springer-Verlag GmbH for their support. I hope you all enjoyed ETAPS 2022.

Finally, a big thanks to Jan, Julia, Dirk, and their local organization team for all their enormous efforts to make ETAPS a fantastic event.

February 2022 Marieke Huisman ETAPS SC Chair ETAPS e.V. President

# Preface

TACAS 2022 was the 28th edition of the International Conference on Tools and Algorithms for the Construction and Analysis of Systems. TACAS 2022 was part of the 25th European Joint Conferences on Theory and Practice of Software (ETAPS 2022), which was held from April 2 to April 7 in Munich, Germany, as well as online due to the COVID-19 pandemic. TACAS is a forum for researchers, developers, and users interested in rigorous tools and algorithms for the construction and analysis of systems. The conference aims to bridge the gaps between different communities with this common interest and to support them in their quest to improve the utility, reliability, flexibility, and efficiency of tools and algorithms for building computer-controlled systems.

There were four submission categories for TACAS 2022:

Papers of categories 1–3 were restricted to 16 pages, and papers of category 4 to six pages.

This year 159 papers were submitted to TACAS, consisting of 112 research papers, five case study papers, 33 regular tool papers, and nine tool demo papers. Authors were allowed to submit up to four papers. Each paper was reviewed by three Program Committee (PC) members, who made use of subreviewers. Similarly to previous years, it was possible to submit an artifact alongside a paper, which was mandatory for regular tool and tool demo papers.

An artifact might consist of a tool, models, proofs, or other data required for validation of the results of the paper. The Artifact Evaluation Committee (AEC) was tasked with reviewing the artifacts based on their documentation, ease of use, and, most importantly, whether the results presented in the corresponding paper could be accurately reproduced. Most of the evaluation was carried out using a standardized virtual machine to ensure consistency of the results, except for those artifacts that had special hardware or software requirements. The evaluation consisted of two rounds. The first round was carried out in parallel with the work of the PC. The judgment of the AEC was communicated to the PC and weighed in their discussion. The second round took place after paper acceptance notifications were sent out; authors of accepted research papers who did not submit an artifact in the first round could submit their artifact at this time. In total, 86 artifacts were submitted (79 in the first round and seven in the second) and evaluated by the AEC regarding their availability, functionality, and/or reusability. Papers with an artifact that was successfully evaluated include one or more badges on the first page, certifying the respective properties.

Selected authors were requested to provide a rebuttal for both papers and artifacts in case a review gave rise to questions. Using the review reports and rebuttals, the Program and the Artifact Evaluation Committees extensively discussed the papers and artifacts and ultimately decided to accept 33 research papers, one case study, 12 tool papers, and four tool demos.

This corresponds to an acceptance rate of 29.46% for research papers and an overall acceptance rate of 31.44%.

Besides the regular conference papers, this two-volume proceedings also contains 16 short papers that describe the participating verification systems and a competition report presenting the results of the 11th SV-COMP, the competition on automatic software verifiers for C and Java programs. These papers were reviewed by a separate Program Committee (PC); each of the papers was assessed by at least three reviewers. A total of 47 verification systems with developers from 11 countries entered the systematic comparative evaluation, including four submissions from industry. Two sessions in the TACAS program were reserved for the presentation of the results: (1) a summary by the competition chair and of the participating tools by the developer teams in the first session, and (2) an open community meeting in the second session.

We would like to thank all the people who helped to make TACAS 2022 successful. First, we would like to thank the authors for submitting their papers to TACAS 2022. The PC members and additional reviewers did a great job in reviewing papers: they contributed informed and detailed reports and engaged in the PC discussions. We also thank the steering committee, and especially its chair, Joost-Pieter Katoen, for his valuable advice. Lastly, we would like to thank the overall organization team of ETAPS 2022.

April 2022 Dana Fisman Grigore Rosu PC Chairs

> Swen Jacobs Andrew Reynolds AEC Chairs, Tools, and Case-study Chairs

> > Dirk Beyer Competition Chair

# Organization

### Program Committee

Saddek Bensalem Verimag, France Nikolaj Bjorner Microsoft, USA

Konstantinos Mamouras Rice University, USA Andrew Reynolds University of Iowa, USA

Parosh Aziz Abdulla Uppsala University, Sweden Luca Aceto Reykjavik University, Iceland Timos Antonopoulos Yale University, USA Dirk Beyer LMU Munich, Germany Jasmin Blanchette Vrije Universiteit Amsterdam, The Netherlands Udi Boker Interdisciplinary Center Herzliya, Israel Hana Chockler King's College London, UK Rance Cleaveland University of Maryland, USA Alessandro Coglio Kestrel Institute, USA Pedro R. D'Argenio Universidad Nacional de Córdoba, Argentina Javier Esparza Technical University of Munich, Germany Bernd Finkbeiner CISPA Helmholtz Center for Information Security, Germany Dana Fisman (Chair) Ben-Gurion University, Israel Martin Fränzle University of Oldenburg, Germany Felipe Gorostiaga IMDEA Software Institute, Spain Susanne Graf Université Joseph Fourier, France Radu Grosu Stony Brook University, USA Arie Gurfinkel University of Waterloo, Canada Klaus Havelund Jet Propulsion Laboratory, USA Holger Hermanns Saarland University, Germany Falk Howar TU Clausthal / IPSSE, Germany Swen Jacobs CISPA Helmholtz Center for Information Security, Germany Ranjit Jhala University of California, San Diego, USA Jan Kretinsky Technical University of Munich, Germany Viktor Kuncak Ecole Polytechnique Fédérale de Lausanne, Switzerland Kim Larsen Aalborg University, Denmark Daniel Neider Max Planck Institute for Software Systems, Germany Dejan Nickovic AIT Austrian Institute of Technology, Austria Corina Pasareanu Carnegie Mellon University, NASA, and KBR, USA Doron Peled Bar Ilan University, Israel Anna Philippou University of Cyprus, Cyprus

# Artifact Evaluation Committee

Fabian Meyer RWTH Aachen University, Germany Stefanie Mohr Technical University of Munich, Germany Malte Mues TU Dortmund, Germany Yuki Nishida Kyoto University, Japan Philip Offtermatt Université de Sherbrooke, Canada Muhammad Osama Eindhoven University of Technology, The Netherlands Jiří Pavela Brno University of Technology, Czech Republic Adrien Pommellet LRDE, France Mathias Preiner Stanford University, USA José Proença CISTER-ISEP and HASLab-INESC TEC, Portugal Tim Quatmann RWTH Aachen University, Germany Etienne Renault LRDE, France Andrew Reynolds (Chair) University of Iowa, USA Mouhammad Sakr University of Luxembourg, Luxembourg Morten Konggaard Schou Aalborg University, Denmark Philipp Schlehuber-Caissier LRDE, France Hans-Jörg Schurr Inria Nancy - Grand Est, France Michael Schwarz Technische Universität München, Germany Joseph Scott University of Waterloo, Canada Ali Shamakhi Tehran Institute for Advanced Studies, Iran Lei Shi University of Pennsylvania, USA Matthew Sotoudeh University of California, Davis, USA Jip Spel RWTH Aachen University, Germany Veronika Šoková Brno University of Technology, Czech Republic

# Program Committee and Jury — SV-COMP

Raveendra Kumar Medicherla

Fatimah Aljaafari University of Manchester, UK Lei Bu Nanjing University, China Thomas Bunk LMU Munich, Germany Marek Chalupa Masaryk University, Czech Republic Priyanka Darke Tata Consultancy Services, India Daniel Dietsch University of Freiburg, Germany Gidon Ernst LMU Munich, Germany Fei He Tsinghua University, China Matthias Heizmann University of Freiburg, Germany Jera Hensel RWTH Aachen University, Germany Falk Howar TU Dortmund, Germany Soha Hussein University of Minnesota, USA Dominik Klumpp University of Freiburg, Germany Henrich Lauko Masaryk University, Czech Republic Will Leeson University of Virginia, USA Xie Li Chinese Academy of Sciences, China Viktor Malík Brno University of Technology, Czech Republic Tata Consultancy Services, India

# Steering Committee

# Additional Reviewers

Abraham, Erika Aguilar, Edgar Akshay, S. Asadi, Sepideh Attard, Duncan Avni, Guy Azeem, Muqsit Bacci, Giorgio Balasubramanian, A. R. Barbanera, Franco Bard, Joachim Basset, Nicolas Bendík, Jaroslav Berani Abdelwahab, Erzana Beutner, Raven Bhandary, Shrajan Biewer, Sebastian

Blicha, Martin Brandstätter, Andreas Bright, Curtis Britikov, Konstantin Brunnbauer, Axel Capretto, Margarita Castiglioni, Valentina Castro, Pablo Ceska, Milan Chadha, Rohit Chalupa, Marek Changshun, Wu Chen, Xiaohong Cruciani, Emilio Dahmen, Sander Dang, Thao Danielsson, Luis Miguel Degiovanni, Renzo Dell'Erba, Daniele Demasi, Ramiro Desharnais, Martin Dierl, Simon Dubslaff, Clemens Egolf, Derek Evangelidis, Alexandros Fedyukovich, Grigory Fiedor, Jan Fitzpatrick, Stephen Fleury, Mathias Frenkel, Hadar Gamboa Guzman, Laura P. Garcia-Contreras, Isabel Gianola, Alessandro Goorden, Martijn Gorostiaga, Felipe Gorrieri, Roberto Grahn, Samuel Grastien, Alban Grover, Kush Grünbacher, Sophie Guha, Shibashis Gutiérrez Brida, Simón Emmanuel Havlena, Vojtěch He, Jie Helfrich, Martin Henkel, Elisabeth Hicks, Michael Hirschkoff, Daniel Hofmann, Jana Hojjat, Hossein Holík, Lukáš Hospodár, Michal Huang, Chao Hyvärinen, Antti Inverso, Omar Itzhaky, Shachar Jaksic, Stefan Jansen, David N. Jin, Xiangyu Jonas, Martin Kanav, Sudeep Karra, Shyam Lal Katsaros, Panagiotis

Kempa, Brian Klauck, Michaela Kreitz, Christoph Kröger, Paul Köhl, Maximilian Alexander König, Barbara Lahijanian, Morteza Larraz, Daniel Le, Nham Lemberger, Thomas Lengal, Ondrej Li, Chunxiao Li, Jianlin Lorber, Florian Lung, David Luppen, Zachary Lybech, Stian Major, Juraj Manganini, Giorgio McCarthy, Eric Mediouni, Braham Lotfi Meggendorfer, Tobias Meira-Goes, Romulo Melcer, Daniel Metzger, Niklas Milovancevic, Dragana Mohr, Stefanie Najib, Muhammad Noetzli, Andres Nouri, Ayoub Offtermatt, Philip Otoni, Rodrigo Paoletti, Nicola Parizek, Pavel Parker, Dave Parys, Paweł Passing, Noemi Perez Dominguez, Ivan Perez, Guillermo Pinna, G. Michele Pous, Damien Priya, Siddharth Putruele, Luciano Pérez, Jorge A. Qu, Meixun Raskin, Mikhail

Rauh, Andreas Reger, Giles Reynouard, Raphaël Riener, Heinz Rogalewicz, Adam Roy, Rajarshi Ruemmer, Philipp Ruijters, Enno Schilling, Christian Schmitt, Frederik Schneider, Tibor Scholl, Christoph Schultz, William Schupp, Stefan Schurr, Hans-Jörg Schwammberger, Maike Shafiei, Nastaran Siber, Julian Sickert, Salomon Singh, Gagandeep Smith, Douglas Somenzi, Fabio

Stewing, Richard Stock, Gregory Su, Yusen Tang, Qiyi Tibo, Alessandro Trefler, Richard Trtík, Marek Turrini, Andrea Vaezipoor, Pashootan van Dijk, Tom Vašíček, Ondřej Vediramana Krishnan, Hari Govind Wang, Wenxi Wendler, Philipp Westfold, Stephen Winter, Stefan Wolovick, Nicolás Yakusheva, Sophia Yang, Pengfei Zeljić, Aleksandar Zhou, Yuhao Zimmermann, Martin

# Contents – Part II

#### Probabilistic Systems

# Contents – Part I

#### Synthesis

David Dill, Wolfgang Grieskamp, Junkil Park, Shaz Qadeer, Meng Xu, and Emma Zhong

#### Grammatical Inference

#### Verification Inference

### Constraint Solving

#### Model Checking and Verification

# Probabilistic Systems

# A Probabilistic Logic for Verifying Continuous-time Markov Chains

Ji Guan<sup>1</sup> and Nengkun Yu<sup>2</sup> ()

<sup>1</sup> State Key Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences, Beijing, China<EMAIL_ADDRESS><sup>2</sup> Centre for Quantum Software and Information, University of Technology Sydney, Sydney, Australia.

<EMAIL_ADDRESS>
Abstract. A continuous-time Markov chain (CTMC) execution is a continuous class of probability distributions over states. This paper proposes a probabilistic linear-time temporal logic, namely continuous-time linear logic (CLL), to reason about the probability distribution execution of CTMCs. We defne the syntax of CLL on the space of probability distributions. The syntax of CLL includes multiphase timed until formulas, and the semantics of CLL allows time reset to study relatively temporal properties. We derive a corresponding model-checking algorithm for CLL formulas. The correctness of the model-checking algorithm depends on Schanuel's conjecture, a central open problem in transcendental number theory. Furthermore, we provide a running example of CTMCs to illustrate our method.

# 1 Introduction

As a popular model of probabilistic continuous-time systems, continuous-time Markov chains (CTMCs) have been extensively studied since Kolmogorov [25]. In the recent 20 years, probabilistic continuous-time model checking receives much attention. Adopting probabilistic computational tree logic (PCTL) [22] to this context with extra multiphase timed until formulas Φ1U <sup>T</sup><sup>1</sup> Φ<sup>2</sup> · · ·U <sup>T</sup><sup>K</sup> ΦK+1, for state formula Φ and time interval T , Aziz et al. proposed continuous stochastic logic (CSL) to specify the branching-time properties of CTMCs and the model-checking problem for CSL is decidable [8]. After that, efcient modelchecking algorithms were developed by transient analysis of CTMCs using uniformization [9] and stratifcation [41] for a restricted version (path formulas are restricted to single until formulas Φ1U <sup>I</sup>Φ2) and a full version of CSL, respectively. These algorithms have been practically implemented in model checkers PRISM [26], MRMC [24] and STORM [18]. Further details can be found in an excellent survey [23].

There are also diferent ways to specify the linear-time properties of CTMCs. Timed automata were frst used to achieve this task [11,13,14,15,19], and then metric temporal logic (MTL) [12] was also considered in this context. Subsequently, the probability of "the system being in state s<sup>0</sup> within fve-time units after having continuously remained in state s1" can be computed. However, some statements cannot be specifed and verifed because of the lack of a probabilistic linear-time temporal logic, for instance "the system being in state s<sup>0</sup> with high probability (≥ 0.9) within fve-time units after having continuously remained in state s<sup>1</sup> with low probability (≤ 0.1)". Furthermore, this probabilistic property cannot be expressed by CSL because CSL cannot express properties that are defned across several state transitions of the same time length in the execution of a CTMC.

In this paper, targeting to express the mentioned probabilistic linear-time properties, we introduce continuous-time linear logic (CLL). In particular, we adopt the viewpoint used in [2] by regarding CTMCs as transformers of probability distributions over states. CLL studies the properties of the probability distribution execution generated by a given initial probability distribution over time. By the fundamental diference between the views of state executions and probability distribution executions of CTMCs, CLL and CSL are incomparable and complementary, as the relation between probabilistic linear-time temporal logic (PLTL) and PCTL in model checking discrete-time Markov chains [2, Section 3.3].

The atomic propositions of CLL are explained on the space of probability distributions over states of CTMCs. We apply the method of symbolic dynamics to the probability distributions of CTMCs. To be specifc, we symbolize the probability value space [0, 1] into a fnite set of intervals I = {I<sup>k</sup> ⊆ [0, 1]} m <sup>k</sup>=1. A probability distribution µ over its set of states S = {s0, s2, . . . , sd−1} is then represented symbolically as a set of symbols

$$\mathbb{S}(\mu) = \{ \langle s, \mathcal{T} \rangle \in S \times \mathcal{J} \,:\, \mu(s) \in \mathcal{T} \}$$

where each symbol ⟨s, I⟩ asserts µ(s) ∈ I, i.e., the probability of state s in distribution µ falls in interval I. For example, ⟨s0, [0.9, 1]⟩ means the system is in state s<sup>0</sup> with a probability in 0.9 to 1. The symbolization idea of distributions has been considered in [2]: choosing a disjoint cover of [0, 1]:

$$\mathcal{T} = \{ [0, p\_1), [p\_1, p\_2), \dots, [p\_n, 1] \}.$$

Here, we remove this restriction and enrich the expressiveness of I . A crucial fact about this symbolization is that the set S × I is fnite. Consequently, the (probability distribution) execution path generated by an initial probability distribution µ induces a sequence of symbols in S × I over time. Therefore, the dynamics of CTMCs can be studied in terms of a (real-time) language over the alphabet S × I , which is the set of atomic propositions of CLL.

Diferent from non-probabilistic linear-time temporal logics — linear-time temporal logic (LTL) and MTL, CLL has two types of formulas: state formulas and path formulas. The state formulas are constructed using propositional connectives. The path formulas are obtained by propositional connectives and a temporal modal operator timed until U T for a bounded time interval T , as in MTL and CSL. The standard next-step temporal operator in LTL is meaningless in continuous-time systems since the time domain (real numbers) is uncountable. As a result, CLL can express the above mentioned probabilistic property "the system is at state s<sup>0</sup> with high probability (≥ 0.9) within 5 time units after having continuously remained at state s<sup>1</sup> with low probability (≤ 0.1)" in a path formula:

$$
\varphi = \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle \dots
$$

In this single until formula, there is a time instant 0 ≤ t ≤ 5 at which state s<sup>1</sup> with low probability transits to state s<sup>0</sup> with high probability. Then we illustrate this on the following timeline.

$$\xleftarrow{\downarrow 0 \overbrace{\langle \ast\_{1}, [0, 0.1] \rangle}^{} \uparrow \langle s\_{0}, [0.9, 1] \rangle}$$

Furthermore, CLL allows multiphase timed until formulas. The semantics of the formulas focuses on relative time intervals, i.e., time can be reset as in timed automata [5,6], while those of CSL [8] are for absolute time intervals. Subsequently, CLL can express not only relatively but also absolutely temporal properties of CTMCs.

We illustrate the signifcant diference between relatively temporal properties and absolutely temporal properties of CTMCs. For instance, "before probability distributions transition φ happening in 3 to 7 time units, the system always stays at state s<sup>0</sup> with a high probability (≥ 0.9)" can be formalized in path formulae

$$
\varphi' = \langle s\_0, [0.9, 1] \rangle U^{[3,7]} (\langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle).
$$

As we can see, there are two time instants, namely t<sup>1</sup> and t2, happening distribution transitions. Time is reset to 0 after the frst distribution transition happens and thus t<sup>2</sup> is relative to t1. More clearly, we depict this on the following timeline.

$$\overbrace{\uparrow0s}^{=3} \underbrace{\downarrow t\_1 \le 7 \qquad \downarrow (t\_2 + t\_1) \le 12}\_{\langle s\_0, [0.9, 1] \rangle \}}^{\downarrow t\_1 \le 7} \langle s\_0, [0.9, 1] \rangle$$

An absolute version is "probability distribution transition φ happens and the system always stays at state s<sup>0</sup> with a high probability (≥ 0.9) in 3 to 7 time units"

$$
\varphi'' = \Box^{[3,7]} \langle s\_0, [0.9, 1] \rangle \wedge \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle \rangle.
$$

We can get a clear timeline representation by simply adding □[3,7]⟨s0, [0.9, 1]⟩ to that of φ. Assume that t < 3,

$$\underbrace{\begin{matrix}\downarrow\\\downarrow\\\downarrow\\\end{matrix}}\_{\langle s\_{1},[0,0.1]\rangle}\underbrace{\begin{matrix}\downarrow\\\langle s\_{0},[0.9,1]\rangle\end{matrix}}\_{\langle\langle s\_{0},[0.9,1]\rangle}\underbrace{\begin{matrix}\downarrow\\\downarrow\\\langle s\_{0},[0.9,1]\rangle\end{matrix}}\_{\langle\langle s\_{0},[0.9,1]\rangle\rangle}\end{matrix}}\_{\langle\langle s\_{0},[0.9,1]\rangle}$$

Time reset enriches the expressiveness of CLL but introduces more difculties to model checking CLL than CSL. We cross this by translating relative time to the absolute one. As a result, we develop an algorithm to model check CTMCs against CLL formulas. More precisely, we reduce the model-checking problem to a reachability problem of absolute time intervals. The reachability problem corresponds to the real root isolation problem of real polynomial-exponential functions (PEFs) over the feld of algebraic numbers, an extensively studied question in recent symbolic and algebraic computation community (e.g. [1,20,28]). By developing a state-of-the-art real root isolation algorithm, we resolve the latter problem under the assumption of the validity of Schanuel's conjecture, a central open question in transcendental number theory [27]. This conjecture has also been the footstone of the correctness of many recent model-checking algorithms, including the decidability of continuous-time Markov decision processes [30], the synthesizing inductive invariants for continuous linear dynamical systems [4], the termination analysis for probabilistic programs with delays [39], and reachability analysis for dynamical systems [20].

In summary, the main contributions of this paper are as follows.

Organization of this paper. In the next section, we give the mathematical preliminaries used in this paper. In Section 3, we recall the view of CTMCs as distribution transformers. After that, the symbolic dynamics of CTMCs are introduced by symbolizing distributions over states of CTMCs in Section 4. In the subsequent section, we present our continuous-time probabilistic temporal logic CLL. In Section 6, we develop an algorithm to solve the CLL model checking problem. A case study and related works are shown in Sections 7 and 8, respectively. We summarize our results and point out future research directions in the fnal section.

#### 2 Preliminaries

For the convenience of the readers, we review basic defnitions and notations of number theory, particularly Schanuel's conjecture.

Throughout this paper, we write C, R, Q and A for the felds of all complex, real, rational and algebraic numbers, respectively. In addition, Z denotes the set of all integer numbers. For F ∈ {C, R, Q, Z, A}, we use F[t] and F <sup>n</sup>×<sup>m</sup> to denote the set of polynomials in t with coefcients in F and n-by-m matrices with every entry in F, respectively. Furthermore, for F ∈ {R, Q, Z}, we use F <sup>+</sup> to denote the set of positive elements (including 0) of F.

A bounded (time) interval T is a subset of R <sup>+</sup>, which may be open, half-open or closed with one of the following forms:

$$[t\_1, t\_2], [t\_1, t\_2), (t\_1, t\_2], (t\_1, t\_2),$$

where t1, t<sup>2</sup> ∈ R <sup>+</sup> and t<sup>2</sup> ≥ t<sup>1</sup> (t<sup>1</sup> = t<sup>2</sup> is only allowed in the case of [t1, t2]). Here, t<sup>1</sup> and t<sup>2</sup> are called the left and right endpoints of T , respectively. Conveniently, we use inf T and sup T to denote t<sup>1</sup> and t2, respectively. In this paper, we only consider bounded intervals.

For reasoning about the temporal properties, we further defne the addition and subtraction of (time) intervals. The expression T + t or t + T , for t ∈ R +, denotes the interval {t + t ′ : t ′ ∈ T }. Similarly, T − t stands for the interval {−t + t ′ : t ′ ∈ T } if t ≤ inf T . Furthermore, for two intervals T<sup>1</sup> and T2,

$$
\mathcal{T}\_1 + \mathcal{T}\_2 = \{ t \in (t' + \mathcal{T}\_2) : t' \in \mathcal{T}\_1 \} = \{ t\_1 + t\_2 : t\_1 \in \mathcal{T}\_1 \text{ and } t\_2 \in \mathcal{T}\_2 \}.
$$

Two intervals T<sup>1</sup> and T<sup>2</sup> are disjoint if their intersection is an empty set, i.e., T<sup>1</sup> ∩ T<sup>2</sup> = ∅. Let us see some concrete examples: 1 + (2, 3) = (3, 4), (2, 3) − 1 = (1, 2), (2, 3) + [3, 4] = (5, 7) and (2, 3), [3, 4] are disjoint. It is obvious that all calculations of time intervals in the above are easy to be computed.

An algebraic number is a complex number that is a root of a non-zero polynomial in one variable with rational coefcients (or equivalent to integer coefcients, by eliminating denominators). An algebraic number α is represented by (P,(a, b), ε) where P is the minimal polynomial of α, a, b ∈ Q and a + bi is an approximation of α such that |α−(a+bi)| < ε and α is the only root of P in the open ball B(a + bi, ε). The minimal polynomial of α is the polynomial with the smallest degree in Q[t] such that α is a root of the polynomial and the coefcient of the highest-degree term is 1. Any root of f(t) ∈ A[t] is algebraic. Moreover, given the representations of a, b ∈ A, the representations of a ± b, <sup>a</sup> b and a · b can be computed in polynomial time, so does the equality checking [17].

Furthermore, a complex number is called transcendental if it is not an algebraic number. In general, it is challenging to verify relationships between transcendental numbers [33]. On the other hand, one can use the Lindemann-Weierstrass theorem to compare some transcendental numbers. The transcendence of e and π are direct corollaries of this theorem.

Theorem 1 (Lindemann-Weierstrass theorem). Let η1, · · · , η<sup>n</sup> be pairwise distinct algebraic complex numbers. Then P k λke <sup>η</sup><sup>k</sup> ̸= 0 for non-zero algebraic numbers λ1, · · · , λn.

The following concepts are introduced to study the general relation between transcendental numbers.

Defnition 1 (Algebraic independence). A set of complex numbers S = {a1, · · · , an} is algebraically independent over Q if the elements of S do not satisfy any nontrivial (non-constant) polynomial equation with coefcients in Q.

By the above defnition, for any transcendental number u, {u} is algebraically independent over Q, while {a} for any algebraic number a ∈ A is not. Thus, a set of complex numbers that is algebraically independent over Q must consist of transcendental numbers. {π, e<sup>π</sup> √ <sup>n</sup>} is also algebraically independent over Q for any positive integer n [31]. Checking the algebraic independence is challenging. For example, it is still widely open whether {e, π} is algebraically independent over Q.

Defnition 2 (Extension feld). Given two felds E ⊆ F, F is an extension feld of E, denoted by F/E, if the operations of E are those of F restricted to E.

For example, under the usual notions of addition and multiplication, the feld of complex numbers is an extension feld of real numbers.

Defnition 3 (Transcendence degree). Let L be an extension feld of Q, the transcendence degree of L over Q is defned as the largest cardinality of an algebraically independent subset of L over Q.

For instance, let Q(e)/Q = {a + be | a, b ∈ Q} and Q( √ 2)/Q = {a + b √ 2 | a, b ∈ Q} be two extension felds of Q. Then the transcendence degree of them are 1 and 0, respectively, by noting that <sup>e</sup> is a transcendental number and <sup>√</sup> 2 is an algebraic number.

Now, Schanuel's conjecture is ready to be presented.

Conjecture 1 (Schanuel's conjecture). Given any complex numbers z1, · · · , z<sup>n</sup> that are linearly independent over Q, the extension feld Q(z1, ..., zn, ez<sup>1</sup> , ..., ez<sup>n</sup> ) has transcendence degree of at least n over Q.

Stephen Schanuel proposed this conjecture during a course given by Serge Lang at Columbia in the 1960s [27]. Schanuel's conjecture concerns the transcendence degree of certain feld extensions of the rational numbers. The conjecture, if proven, would generalize the most well-known results in transcendental number theory signifcantly [29,37]. For example, the algebraical independence of {e, π} would simply follow by setting z<sup>1</sup> = 1 and z<sup>2</sup> = πi, and using Euler's identity e πi + 1 = 0.

# 3 Continuous-time Markov Chains as Distributions Transformers

We begin with the defnition of continuous-time Markov chains (CTMCs). A CTMC is a Markovian (memoryless) stochastic process that takes values on a fnite state set S (|S| = d < ∞) and evolves in continuous-time t ∈ R <sup>+</sup>. Formally,

Defnition 4. A CTMC is a pair M = (S, Q), where S (|S| = d) is a fnite state set and Q ∈ Qd×<sup>d</sup> is a transition rate matrix.

A transition rate matrix Q is a matrix whose of-diagonal entries {Qi,j}i̸=<sup>j</sup> are nonnegative rational numbers, representing the transition rate from state i to state j, while the diagonal entries {Qj,j} are constrained to be − P <sup>j</sup≯=<sup>i</sup> Qi,j for all 1 ≤ j ≤ d. Consequently, the column summations of Q are all zero.

The evolution of a CTMC can be regarded as a distribution transformer. Given initial distribution µ ∈ Qd×<sup>1</sup> ∈ D(S), the distribution at time t ∈ R <sup>+</sup> is:

$$
\mu\_t = e^{Qt} \mu,
$$

where D(S) is denoted as the set of all probability distributions over S. We call D(S) the probability distribution space of CTMCs. An execution path of CTMCs is a continuous function indexed by initial distribution µ ∈ D(S):

$$
\sigma\_{\mu} \colon \mathbb{R}^+ \to \mathcal{D}(\mathcal{S}), \qquad \sigma\_{\mu}(t) = e^{Qt}\mu. \tag{1}
$$

Example 1. We recall the illustrating example of CTMC M = (S, Q) in [8, Figure 1] as the running example in our work. In particular, M is a 5-dimensional CTMC with initial distribution µ, where S = {s0, s1, s2, s3, s4} and

#### 4 Symbolic Dynamics of CTMCs

In this section, we introduce symbolic dynamics to characterize the properties of the probability distribution space of CTMCs.

First, we fx a fnite set of intervals I = {I<sup>k</sup> ⊆ [0, 1]}k∈K, where the endpoints of each I<sup>k</sup> are rational numbers. With the states S = {s0, s1, · · · , sd−1}, we defne the symbolization of distributions as a function:

$$\mathbb{S}: \mathcal{D}(\mathcal{S}) \to 2^{S \times \mathcal{J}} \qquad \mathbb{S}(\mu) = \{ \langle s, \mathcal{Z} \rangle \in S \times \mathcal{J} \,:\, \mu(s) \in \mathcal{Z} \},\tag{2}$$

where × denotes the Cartesian product, and 2<sup>S</sup>×<sup>I</sup> is the power set of S × I . ⟨s, I⟩ ∈ S(µ) asserts that the probability of state s in distribution µ is in the interval I. The symbolization of distributions is a generalization of the discretization of distributions with Ik∩I<sup>m</sup> = ∅ for all k ̸= m which was studied in [2]. This generalization increases the expressiveness of our continuous linear-time logic introduced in the next section. Now, we can represent any given probability distribution by fnite symbols from S × I . For example, suppose

$$\mathcal{A}' = \{ [0, 0.1], (0.1, 0.9), [0.9, 1], [1, 1], [0.4, 0.4] \}, \tag{3}$$

and then the initial distribution µ in Example 1 is symbolized as

$$\begin{aligned} \mathbb{S}(\mu) &= \{ \langle s\_0, [0, 0.1] \rangle, \langle s\_1, (0.1, 0.9) \rangle, \langle s\_2, (0.1, 0.9) \rangle, \\ &\quad \langle s\_3, (0.1, 0.9) \rangle, \langle s\_3, [0.4, 0.4] \rangle, \langle s\_4, [0, 0.1] \rangle \}. \end{aligned} \tag{4}$$

As we can see from the above example, the symbolization of distributions on states considers the exact probabilities (singleton intervals) of the states and the range of their possibilities.

Next, we introduce the symbolization to CTMCs,

Defnition 5. A symbolized CTMC is a tuple SM = (S, Q, I ), where M = (S, Q) is a CTMC and I is a fnite set of intervals in [0, 1].

As we can see, the set of intervals is picked depending on CTMCs. Then, we extend this symbolization to the path σµ:

$$\mathbb{S}\diamond \sigma\_{\mu} : \mathbb{R}^{+} \to 2^{S \times \mathcal{J}} . \tag{5}$$

Defnition 6. Given a symbolized CTMC SM = (S, Q, I ), S◦σ<sup>µ</sup> is a symbolic execution path of M = (S, Q).

Given a symbolized CTMC SM = (S, Q, I ), the path σ<sup>µ</sup> of CTMC M = (S, Q) over real numbers R <sup>+</sup> generated by probability distribution µ induces a symbolic execution path S ◦ σ<sup>µ</sup> over fnite symbols S × I . Subsequently, the dynamics of CTMCs can be studied in terms of a language over S × I . In other words, we can study the temporal properties of CTMCs in the context of symbolized CTMCs.

# 5 Continuous Linear-time Logic

In this section, we introduce continuous linear-time logic (CLL), a probabilistic linear-time temporal logic, to specify the temporal properties of a symbolized CTMC SM = (S, Q, I ).

CLL has two types of formulas: state formulas and path formulas. The state formulas are constructed using propositional connectives. The path formulas are obtained by propositional connectives and a temporal modal operator timed until U T for a bounded time interval T , as in MTL and CSL. Furthermore, multiphase timed until formulas Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> . . . U<sup>T</sup><sup>n</sup> Φ<sup>n</sup> are allowed to enrich the expressiveness of CLL. More importantly, time reset is involved in these multiphase formulas. Thus absolutely and relatively temporal properties of CTMCs can be studied.

Defnition 7. The state formulas of CLL are described according to the following syntax:

$$\Phi := \mathbf{true} \mid a \in AP \mid \neg \Phi \mid \Phi\_1 \land \Phi\_2$$

where AP denotes S × I as the set of atomic propositions.

The path formulas of CLL are constructed by the following syntax:

φ := true | Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> . . . U<sup>T</sup><sup>n</sup> Φ<sup>n</sup> | ¬φ | φ<sup>1</sup> ∧ φ<sup>2</sup>

where n ∈ Z <sup>+</sup> is a positive integer, for all 0 ≤ k ≤ n, Φ<sup>k</sup> is a state formula, and Tk's are time intervals with the endpoints in Q<sup>+</sup>, i.e., each T<sup>k</sup> is one of the following forms:

(a, b), [a, b],(a, b], [a, b) ∀a, b ∈ Q +.

The semantics of CLL state formulas is defned on the set D(S) of probability distributions over S with the symbolized function S in Eq.(2) of Section 4.

The semantics of CLL path formulas is defned on execution paths {σµ}µ∈D(S) of CTMC M = (S, Q).

Not surprisingly, other Boolean connectives are derived in the standard way, i.e., false = ¬true, Φ<sup>1</sup> ∨ Φ<sup>2</sup> = ¬(¬Φ<sup>1</sup> ∧ ¬Φ2) and Φ<sup>1</sup> → Φ<sup>2</sup> = ¬Φ<sup>1</sup> ∨ Φ2, and the path formula φ follows the same way. Furthermore, we generalize temporal operators ♢ ("eventually") and □ ("always") of discrete-time systems into their timed variant ♢ <sup>T</sup> and □<sup>T</sup> , respectively, in the following:

$$
\Diamond^{\top} \Phi = \mathbf{tr} \mathbf{u} \mathbf{e} U^{\top} \Phi \qquad \Box^{\top} \Phi = \neg \Diamond^{\top} \neg \Phi.
$$

For n = 1 in multiphase timed until formulas, the until operator U T<sup>1</sup> is a timed variant of the until operator of LTL; the path formula Φ0U <sup>T</sup><sup>1</sup> Φ<sup>1</sup> asserts that Φ<sup>1</sup> is satisfed at some time instant in the interval T<sup>1</sup> and that at all preceding time instants in T1, Φ<sup>0</sup> holds. For example,

$$
\varphi = \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle,
$$

as mentioned in introduction section.

For general n, the CLL path formula Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> . . . U<sup>T</sup><sup>n</sup> Φ<sup>n</sup> is explained over the induction on n. We frst mention that U T is right-associative, e.g., Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> stands for Φ0U <sup>T</sup><sup>1</sup> (Φ1U <sup>T</sup><sup>2</sup> Φ2). This makes time reset, i.e., T<sup>1</sup> and T<sup>2</sup> do not have to be disjoint, and the starting time point of T<sup>2</sup> is based on some time instant in T1. Recall the multiphase timed until formula in introduction section and this formula expresses a relative time property:

$$
\varphi' = \langle s\_0, [0.9, 1] \rangle U^{[3,7]} (\langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle),
$$

which is diferent to the following CLL path formula representing an absolutely temporal property of CTMCs:

$$
\varphi'' = \Box^{[3,7]} \langle s\_0, [0.9, 1] \rangle \wedge \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle \rangle.
$$

As an example, we clarify the semantics of CLL by comparing the above two path formulas in general forms:

$$
\Phi\_0 U^{\mathcal{T}\_1} \Phi\_1 U^{\mathcal{T}\_2} \Phi\_2 \quad \text{and} \quad \Phi\_0 U^{\mathcal{T}\_1} \Phi\_1 \wedge \Phi\_1 U^{\mathcal{T}\_2} \Phi\_2 .
$$

(1) σ<sup>µ</sup> |= Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> asserts that there are time instants t<sup>1</sup> ∈ T1, t<sup>2</sup> ∈ T<sup>2</sup> such that µt1+t<sup>2</sup> |= Φ<sup>2</sup> and for any t ′ <sup>1</sup> ∈ T<sup>1</sup> ∩ [0, t1) and t ′ <sup>2</sup> ∈ T<sup>2</sup> ∩ [0, t2), µt ′ 1 |= Φ<sup>0</sup> and µt1+<sup>t</sup> ′ 2 |= Φ1, where µ<sup>t</sup> = e Qtµ ∀t ∈ R <sup>+</sup>. This is more clear in the following timeline.

$$\underbrace{\overbrace{\uparrow\text{time }0}^{\text{inf }\mathscr{T}\_{1}}}\_{\uparrow\text{ time }0} \underbrace{\overbrace{\uparrow t\_{1} \le \sup\,\mathscr{T}\_{1}}^{\text{inf }\mathscr{T}\_{2}}}\_{\Phi\_{0}} \downarrow \Phi\_{2}$$

(2)  $\sigma\_{\mu} = \Phi\_0 U^{\mathcal{T}\_1} \Phi\_1 \wedge \Phi\_1 U^{\mathcal{T}\_2} \Phi\_2$  asserts that there are time instants  $t\_1 \in \mathcal{T}\_1, t\_2 \in \mathcal{T}\_2$  such that  $\mu\_{t\_1} = \Phi\_1$  and  $\mu\_{t\_2} = \Phi\_2$ , and for any  $t\_1' \in \mathcal{T}\_1 \cap [0, t\_1)$  and  $t\_2' \in \mathcal{T}\_2 \cap [0, t\_2)$ ,  $\mu\_{t\_1'} = \Phi\_0$  and  $\mu\_{t\_2'} = \Phi\_1$ , where  $\mu\_t = e^{Qt} \mu$   $\forall t \in \mathbb{R}^+$ .

Before solving the model-checking problem of CTMCs against CLL formulas in the next section, we shall frst discuss what can be specifed in our logic CLL.

Given a CTMC (S, Q), CLL path formula ♢ [0,1000]⟨s, [1, 1]⟩ expresses a liveness property that state s ∈ S is eventually reached with probability one before time instant 1000. In terms of safety properties, formula □[100,1000]⟨s, [0, 0]⟩ represents that state s ∈ S is never reached (reached with probability zero) between time instants 100 and 1000. Furthermore, setting the intervals nontrivial (neither [0, 0] or [1, 1]), liveness and safety properties can be asserted with probabilities, such as ♢ [0,1000]⟨s, [0.5, 1]⟩ and □[100,1000]⟨s, [0, 0.5]⟩. For multiphase timed until formula ⟨s, [0.7, 1]⟩U [2,3]⟨s, [0.7, 1]⟩. . . U[2,3]⟨s, [0.7, 1]⟩, where the number of U [2,3] is 100, asserts that the probability of state s is beyond 0.7 in every time instant 2 to 3, and this happens at least 100 times.

Next, we can classify members of I as representing "low" and "high" probabilities. For example, if I contains 3 intervals {[0, 0.1],(0.1, 0.9), [0.9, 1]}, we can declare the frst interval as "low" and the last interval as "high". In this case □[10,1000)(⟨s0, [0, 0.1]⟩ → ⟨s1, [0.9, 1]⟩) says that, in time interval [10, 1000), whenever the probability of state s<sup>0</sup> is low, the probability of state s<sup>1</sup> will be high.

#### 6 CLL Model Checking

In this section, we provide an algorithm to model check CTMCs against CLL formulas, i.e., the following CLL model-checking problem — Problem 1 is decidable.

Problem 1 (CLL Model-checking Problem). Given a symbolized CTMC SM = (S, Q, I ) with an initial distribution µ and a CLL path formula φ on AP = S ×I , the goal is to decide whether σ<sup>µ</sup> |= φ, where σµ(t) = e Qtµ is an execution path defned in Eq.(1).

In particular, we show that

Theorem 2. Under the condition that Schanuel's conjecture holds, the CLL model-checking problem in Problem 1 is decidable.

In the following, we prove the above theorem from checking basic formulas — atomic propositions to the most complex one — nontrivial multiphase timed until formulas. For readability, we put the proofs of all results in Appendix A of the extended version [21] of this paper.

We start with the simplest case of atomic proposition ⟨s, I⟩. By the semantics of CLL, µ<sup>t</sup> |= ⟨s, I⟩ if and only if µ<sup>t</sup> = e Qtµ(s) ∈ I. To check this, we frst observe that the execution path e Qtµ of CTMCs is a system of polynomial exponential functions (PEFs).

Defnition 8. A function f : R → R is a polynomial-exponential function (PEF) if f has the following form:

$$f(t) = \sum\_{k=0}^{K} f\_k(t)e^{\lambda\_k t} \tag{6}$$

where for all 0 ≤ k ≤ K < ∞, fk(t) ∈ F1[t], fk(t) ̸= 0, λ<sup>k</sup> ∈ F<sup>2</sup> and F1, F<sup>2</sup> are felds. Without loss of generality, we assume that λk's are distinct.

Generally, for a PEF f(t) with the range in complex numbers C, g(t) = f(t) + f ∗ (t) is a PEF with the range in real numbers R, where f ∗ (t) is the complex conjugate of f(t). The factor t is omitted whenever convenient, i.e., f = f(t). t is called a root of a function f if f(t) = 0. PEFs often appear in transcendental number theory as auxiliary functions in the proofs involving the exponential function [10].

Lemma 1. Given a CTMC M = (S, Q) with S = {s0, . . . , sd−1}, Q ∈ Qd×<sup>d</sup> , and an initial distribution µ ∈ Qd×<sup>1</sup> , for any 0 ≤ i ≤ d − 1 , e Qtµ(si), the i-th entry of e Qtµ, can be expressed as a PEF f : R <sup>+</sup> → [0, 1] as in Eq.(6) with F<sup>1</sup> = F<sup>2</sup> = A.

By the above lemma, for a given t in some bounded time interval T (to be specifc in the latter discussion), e Qtµ(s) ∈ I is determined by the algebraic structure of PEF g(t) = e Qtµ(s) in T . That is all maximum intervals Tmax ⊆ T such that g(t) ∈ I for all t ∈ Tmax, where interval Tmax ̸= ∅ is called maximum for g(t) ∈ I if no sub-intervals T ′ ⊊ Tmax such that the property holds, i.e., g(t) ∈ I for all t ∈ T ′ . Then e Qtµ(s) ∈ I if and only if t ∈ Tmax for some maximum interval Tmax. So, we aim to compute the set T of all maximum intervals. By the continuity of PEF g(t), this can be done by identifying a real root isolation of the following PEF f(t) in T : f(t) = (g(t) − inf I)(g(t) − sup I).

A (real) root isolation of function f(t) in interval T is a set of mutually disjoint intervals, denoted by Iso(f)<sup>T</sup> = {(a<sup>j</sup> , b<sup>j</sup> ) ⊆ T } for a<sup>j</sup> , b<sup>j</sup> ∈ Q such that

– for any j, there is one and only one root of f(t) in (a<sup>j</sup> , b<sup>j</sup> );

– for any root t <sup>∗</sup> of f(t), t <sup>∗</sup> ∈ (a<sup>j</sup> , b<sup>j</sup> ) for some j.

Furthermore, if f has no any root in T , then Iso(f)<sup>T</sup> = ∅.

Although there are infnite kinds of real root isolations of f(t) in T , the number of isolation intervals equals to the number of distinct roots of f(t) in T .

Finding real root isolations of PEFs is a long-standing problem and can be at least backtracked to Ritt's paper [34] in 1929. Some following results were obtained since the last century (e.g. [7,38]). This problem is essential in the reachability analysis of dynamical systems, one active feld of symbolic and algebraic computation. In the case of F<sup>1</sup> = Q and F<sup>2</sup> = N <sup>+</sup> in [1], an algorithm named ISOL was proposed to isolate all real roots of f(t). Later, this algorithm has been extended to the case of F<sup>1</sup> = Q and F<sup>2</sup> = R [20]. A variant of the problem has also been studied in [28]. The correctness of these algorithms is based on Schanuel's conjecture. Other works are using Schanuel's conjecture to do the root isolation of other functions, such as exp-log functions [35] and tame elementary functions [36].

By Lemma 1, we pursue this problem in the context of CTMCs. The distinct feature of solving real root isolations of PEFs in our paper is to deal with complex numbers C, more specifcally algebraic numbers A, i.e., F<sup>1</sup> = F<sup>2</sup> = A. At the same time, to the best of our knowledge, all the previous works can only handle the case over R. Here, we develop a state-of-the-art real root isolation algorithm for PEFs over algebraic numbers. Thus from now on, we always assume that PEFs are over A, i.e., F<sup>1</sup> = F<sup>2</sup> = A in Eq.(6). In this case, it is worth noting that whether a PEF has a root in a given interval, T ⊆ R <sup>+</sup> is decidable subject to Schanuel's Conjecture if T is bounded [16], which falls in the situation we consider in this paper.

Theorem 3 ([16]). Under the condition that Schanuel's conjecture holds, there is an algorithm to check whether a PEF f(t) has a root in interval T , i.e., whether Iso(f)<sup>T</sup> = ∅.

In this paper, we extend the above checking Iso(f)<sup>T</sup> = ∅ to computing Iso(f)<sup>T</sup> of PEF f(t).

Theorem 4. Under the condition that Schanuel's conjecture holds, there is an algorithm to fnd real root isolation Iso(f)<sup>T</sup> for any PEF f(t) and interval T . Furthermore, the number of real roots is fnite, i.e., |Iso(f)<sup>T</sup> | < ∞.

We can compute the set T of all maximum intervals with the above theorem to check atomic propositions. Furthermore, we can compare the values of any real roots of PEFs, which is important in model checking general multiphase timed until formulas at the end of this section.

Lemma 2. Let f1(t) and f2(t) be two PEFs with the domains in T<sup>1</sup> and T2, and t<sup>1</sup> ∈ T<sup>1</sup> and t<sup>2</sup> ∈ T<sup>2</sup> are roots of them, respectively. Under the condition that Schanuel's conjecture holds, there is an efcient way to check whether or not t<sup>1</sup> − t<sup>2</sup> < g for any given rational number g ∈ Q.

For model checking general state formula Φ, we can also use real root isolation of some PEF to obtain the set of all maximum intervals Tmax such that µ<sup>t</sup> |= Φ for all t ∈ Tmax. The reason is that Φ admits conjunctive normal form consisting of atomic propositions. See the proof of the following lemma in Appendix A of the extended version [21] of this paper for the details.

Lemma 3. Under the condition that Schanuel's conjecture holds, given a time interval T , the set T of all maximum intervals in T satisfying µ<sup>t</sup> |= Φ can be computed, where Φ is a state formula of CLL. Furthermore, the number of all intervals in T is fnite; the left and right endpoints of each interval in T are roots of PEFs.

At last, we characterize the multiphase timed until formulas by the reachability analysis of time intervals (instants).

Lemma 4. σ<sup>µ</sup> |= Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> · · ·U <sup>T</sup><sup>n</sup> Φ<sup>n</sup> if and only if there exist time intervals {I<sup>k</sup> ⊆ R <sup>+</sup>} n <sup>k</sup>=0 with I<sup>0</sup> = [0, 0] such that

By the above lemma, the problem of checking multiphase timed until formulas is reduced to verify the existence of a sequence of time intervals.

Now we can show the proof of Theorem 2.

Proof. Recall that the nontrivial step is to model check multiphase timed until formula Φ0U <sup>T</sup><sup>1</sup> Φ1U <sup>T</sup><sup>2</sup> Φ<sup>2</sup> · · ·U <sup>T</sup><sup>n</sup> Φn, where {Tj} n <sup>j</sup>=1 is a set of bounded rational intervals in R <sup>+</sup>, and for 0 ≤ k ≤ n + 1, Φ<sup>k</sup> is a state formula.

By Lemma 4, for model checking the above formula, we only need to check the existence of time intervals {Ik} n <sup>k</sup>=0 illustrated in the lemma. The following procedure can construct such a set of intervals if it exists:

$$\mathcal{J}\_k = \{ \mathcal{T} \cap (\mathcal{T}' + \mathcal{T}\_k) : \mathcal{T} \in \mathcal{J}\_k \text{ and } \mathcal{T}' \in \mathcal{J}\_{k-1} \}. \tag{7}$$

The above updates can be fnished by Lemma 2. If I<sup>k</sup> = ∅, then the formula is not satisfed;

– (4) Updating In: for each I ∈ In, we replace I with [s − ε, s) for some constant ε > 0 if there is an s ∈ I with s − ε ∈ I such that µ<sup>s</sup> |= Φ<sup>n</sup> where µ<sup>s</sup> = e Qsµ; Otherwise, remove this element from In. Again, this can be done by Lemma 3. If I<sup>n</sup> = ∅, then the formula is not satisfed;

– (5) Finally, let k from n − 1 to 1, updating Ik:

$$\mathcal{A}\_k^\prime = \{ [s-\inf \mathcal{T}\_k, s-\inf \mathcal{T}\_k] : [s-\varepsilon, s) \in \mathcal{A}\_{k+1}^\prime \}.$$

Thus after the above procedure, we have non-empty sets {Ik} n <sup>k</sup>=0 with the following properties.

Therefore, we can get a set of intervals {Ik} n <sup>k</sup>=0 satisfying the two conditions in Lemma 4 if it exists. On the other hand, it is easy to check that all such {Ik} n <sup>k</sup>=0 must be in {Ik} n <sup>k</sup>=0, i.e., for each k, I<sup>k</sup> ⊆ I for some I ∈ Ik. This ensures the correctness of the above procedure.

By the above constructive analysis, we give an algorithm for model checking CTMCs against CLL formulas. Focusing on the decidability problem, we do not provide the pseudocode of the algorithm. Alternatively, we implement a numerical experiment to illustrate the checking procedure in the next section.

#### 7 Numerical Implementation

In this section, we implement a case study of checking CTMCs against CLL formulas. Here, we consider a symbolized CTMC SM = (S, Q, I ), where M = (S, Q) is the CTMC in Example 1 and fnite set I is the one considered in Eq.(3). We check the properties of M given by the following two CLL path formulas mentioned in the introduction for diferent initial distributions.

$$\begin{aligned} \varphi &= \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle . \\ \varphi' &= \langle s\_0, [0.9, 1] \rangle U^{[3,7]} \langle s\_1, [0, 0.1] \rangle U^{[0,5]} \langle s\_0, [0.9, 1] \rangle . \end{aligned}$$

By Jordan decomposition, we have Q = SJS<sup>−</sup><sup>1</sup> where

$$S = \begin{pmatrix} 0 & -6 \ 0 \ 0 \ 0 \\ 0 & 2 \ 0 \ 0 \ 1 \\ -7 - 3 \ 0 \ 0 \ 0 \\ 3 & 3 \ 0 \ 1 \ 0 \\ 4 & 4 & 1 \ 0 \ 0 \end{pmatrix} \qquad J = \begin{pmatrix} -7 & 0 \ 0 \ 0 \ 0 \\ 0 & -3 \ 0 \ 0 \ 0 \\ 0 & 0 \ 0 \ 0 \ 0 \\ 0 & 0 \ 0 \ 0 \ 0 \\ 0 & 0 \ 0 \ 0 \ 0 \end{pmatrix} \qquad S^{-1} = \begin{pmatrix} \frac{1}{14} & 0 - \frac{1}{7} & 0 & 0 \\ -\frac{1}{6} & 0 & 0 & 0 \\ \frac{8}{21} & 0 & \frac{4}{7} & 0 \ 1 \\ \frac{2}{7} & 0 & \frac{3}{7} & 1 \ 0 \\ \frac{1}{3} & 1 & 0 & 0 \ 0 \end{pmatrix}.$$

Then, we consider an initial distribution µ as the same as the one in Example 1. Then we have that the value of e Qtµ is as follows:

$$
\begin{pmatrix} e^{-3t} & 0 & 0 & 0 \\ -\frac{1}{3}(e^{-3t}-1) & 1 & 0 & 0 \\ \frac{1}{2}(e^{-3t}-e^{-7t}) & 0 & e^{-7t} & 0 \\ \frac{3}{14}e^{-7t}-\frac{1}{2}e^{-3t}+\frac{2}{7}\ 0 & -\frac{3}{7}e^{-7t}+\frac{3}{7}\ 1 & 0 \\ \frac{1}{7}e^{-7t}-\frac{2}{3}e^{-3t}+\frac{8}{21}\ 0 & -\frac{4}{7}e^{-7t}+\frac{4}{7}\ 0 & 1 \\ \end{pmatrix}
\begin{pmatrix} 0.1 \\ 0.2 \\ 0.3 \\ 0.4 \\ 0 \end{pmatrix} = \begin{pmatrix} \frac{1}{10}e^{-3t} \\ -\frac{1}{30}e^{-3t}+\frac{7}{30} \\ \frac{1}{20}e^{-3t}+\frac{1}{4}e^{-7t} \\ -\frac{1}{20}e^{-3t}-\frac{3}{28}e^{-7t}+\frac{39}{20} \\ -\frac{1}{15}e^{-3t}-\frac{7}{7}e^{-7t}+\frac{39}{105} \end{pmatrix}.
$$

As we only consider states s<sup>0</sup> and s<sup>1</sup> in formulas φ and φ ′ , we focus on the following PEFs: f0(t) = <sup>1</sup> <sup>10</sup> e <sup>−</sup>3<sup>t</sup> and f1(t) = − 1 <sup>30</sup> e <sup>−</sup>3<sup>t</sup> + 7 30 .

Next, we initialize the model checking procedures introduced in the proof of Theorem 2. First, we compute the set T of all maximum intervals T ⊆ [0, 5] such that e Qtµ |= ⟨s0, [0.9, 1]⟩ for t ∈ T , i.e., f0(t) ∈ [0.9, 1] for t ∈ T . We obtain T = ∅ by the real root isolation algorithm mentioned in Theorem 4, and this indicates that σ<sup>µ</sup> ̸|= φ where σµ(t) = e Qtµ is the path induced by µ and defned in Eq.(1).

To check whether σ<sup>µ</sup> |= φ ′ , we compute the set T of all maximum intervals T ⊆ [0, 12] such that e Qtµ |= ⟨s0, [0.9, 1]⟩ for t ∈ T , i.e., f0(t) ∈ [0.9, 1] for t ∈ T . Again, we obtain T = ∅ by the real root isolation algorithm in Theorem 4. Therefore, σ<sup>µ</sup> ̸|= φ ′ .

In the following, we consider a diferent initial distribution µ<sup>1</sup> as follows:

$$e^{Qt}\mu\_1 = e^{Qt} \begin{pmatrix} 0.9\\0\\0.1\\0\\0 \end{pmatrix} = \begin{pmatrix} \frac{9}{10}e^{-3t} \\ -\frac{3}{10}(e^{-3t}-1) \\\ \frac{9}{20}e^{-3t}-\frac{7}{20}e^{-7t} \\ -\frac{9}{20}e^{-3t}+\frac{3}{20}e^{-7t}+\frac{3}{10} \\ -\frac{3}{5}e^{-3t}+\frac{1}{5}e^{-7t}+\frac{2}{5} \end{pmatrix}.$$

The key PEFs are: g0(t) = <sup>9</sup> <sup>10</sup> e <sup>−</sup>3<sup>t</sup> and g1(t) = − 3 <sup>10</sup> (e <sup>−</sup>3<sup>t</sup> − 1).

Again, we initialize the model checking procedures introduced in the proof of Theorem 2. We frst compute the set T of all maximum intervals T ⊆ [0, 5] such that e Qtµ<sup>1</sup> |= ⟨s1, [0, 0.1]⟩ for t ∈ T , i.e., g1(t) ∈ [0, 0.1] for t ∈ T . This can be done by fnding a real root isolation of the following PEF: g 0 1 (t) = − 3 <sup>10</sup> (e <sup>−</sup>3<sup>t</sup> − 1) − 1 10 .

By implementing the real root isolation algorithm in Theorem 4, we have

$$\text{Iso}(g\_1^0)\_{[0,5]} = \{ (0.13, 0.14) \} \text{ and then } \mathcal{F} = \{ [0, t^\*] \} \text{ for } t^\* \in (0.13, 0.14).$$

Following the same way, we compute T for e Qtµ<sup>1</sup> |= ⟨s0, [0.9, 1]⟩. Then we complete the model checking procedures in the proof of Theorem 2, and we conclude: σ<sup>µ</sup><sup>1</sup> |= φ. By repeating these, the result of the second formula φ ′ is σ<sup>µ</sup><sup>1</sup> ̸|= φ ′ .

#### 8 Related Works

Agrawal et al. [2] introduced probabilistic linear-time temporal logic (PLTL) to reason about discrete-time Markov chains in the context of distribution transformers as we did for CTMCs in this paper. Interestingly, the Skolem Problem can be reduced to the model checking problem for the logic PLTL [3]. The Skolem Problem asks whether a given linear recurrence sequence has a zero term and plays a vital role in the reachability analysis of linear dynamical systems. Unfortunately, the decidability of the problem remains open [32]. Recently, the Continuous Skolem Problem has been proposed with good behavior (the problem is decidable) and forms a fundamental decision problem concerning reachability

in continuous-time linear dynamical systems [16]. Not surprisingly, the Continuous Skolem Problem can be reduced to model-checking CLL. The primary step of verifying CLL formulas is to fnd a real root isolation of a PEF in a given interval. Chonev, Ouaknine and Worrell reformulated the Continuous Skolem Problem in terms of whether a PEF has a root in a given interval, which is decidable subject to Schanuel's conjecture [16]. An algorithm for fnding root isolation can also answer the problem of checking the existence of the roots of a PEF. However, the reverse does not work in general. Therefore, the decidability of the Continuous Skolem Problem cannot be applied to establish that of our CLL model checking.

Remark 1. By adopting the method in this paper, we established the decidability of model checking quantum CTMCs against signal temporal logic [40]. Again, we need Schanuel's conjecture to guarantee the correctness. A Lindblad's master equation governs a quantum CTMC and a more general real-time probabilistic Markov model than a CTMC, i.e., a CTMC is an instance of quantum CTMCs. We converted the evolution of Lindblad's master equation into a distribution transformer that preserves the laws of quantum mechanics. We reduced the model-checking problem of quantum CTMCs to the real root isolation problem, which we considered in this paper, and thus our method could be applied to it.

# 9 Conclusion

This paper revisited the study of temporal properties of fnite-state CTMCs by symbolizing the probability value space [0, 1] into a fnite set of intervals. To specify relatively and absolutely temporal properties, we propose a probabilistic logic for CTMCs, namely continuous linear-time logic (CLL). We have considered the model checking problem in this setting. Our main result is that a state-of-theart real root isolation algorithm over the feld of algebraic numbers was proposed to establish the decidability of the model checking problem under the condition that Schanuel's conjecture holds.

This paper aims to show decidability in as simple a fashion as possible without paying much attention to complexity issues. Faster algorithms on our current constructions would signifcantly improve from a practical standpoint.

# Acknowledgments

We want to thank Professor Joost-Pieter Katoen for his invaluable feedback and for pointing out the references [14,15,30]. This work is supported by the National Key R&D Program of China (Grant No: 2018YFA0306701), the National Natural Science Foundation of China (Grant No: 61832015), ARC Discovery Program (#DP210102449) and ARC DECRA (#DE180100156).

#### References

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# Under-Approximating Expected Total Rewards in POMDPs?

Alexander Bork1() , Joost-Pieter Katoen<sup>1</sup> , and Tim Quatmann<sup>1</sup>

> RWTH Aachen University, Aachen, Germany<EMAIL_ADDRESS>
Abstract We consider the problem: is the optimal expected total reward to reach a goal state in a partially observable Markov decision process (POMDP) below a given threshold? We tackle this—generally undecidable—problem by computing under-approximations on these total expected rewards. This is done by abstracting finite unfoldings of the infinite belief MDP of the POMDP. The key issue is to find a suitable under-approximation of the value function. We provide two techniques: a simple (cut-off) technique that uses a good policy on the POMDP, and a more advanced technique (belief clipping) that uses minimal shifts of probabilities between beliefs. We use mixed-integer linear programming (MILP) to find such minimal probability shifts and experimentally show that our techniques scale quite well while providing tight lower bounds on the expected total reward.

# 1 Introduction

The relevance of POMDPs. Partially observable Markov decision processes (POM-DPs) originated in operations research and nowadays are a pivotal model for planning in AI [40]. They inherit all features of classical MDPs: each state has a set of discrete probability distributions over the states and rewards are earned when taking transitions. However, states are not fully observable. Intuitively, certain aspects of the states can be identified, such as a state's colour, but states themselves cannot be observed. This partial observability reflects, for example, a robot's view of its environment while only having the limited perspective of its sensors at its disposal. The main goal is to obtain a policy—a plan how to resolve the non-determinism in the model—for a given objective. The key problem here is that POMDP policies must base their decisions only on the observable aspects (e.g. colours) of states. This stands in contrast to policies for MDPs which can make decisions dependent on the entire history of full state information.

Analysing POMDPs. Typical POMDP planning problems consider either finitehorizon objectives or infinite-horizon objectives under discounting. Finite-horizon objectives focus on reaching a certain goal state (such as "the robot has collected all items") within a given number of steps. For infinite horizons, no step bound

<sup>?</sup> This work is funded by the DFG RTG 2236 "UnRAVeL".

is provided and typically rewards along a run are weighted by a discounting factor that indicates how much immediate rewards are favoured over more distant ones. Existing techniques to treat these objectives include variations of value iteration [46,36,20,18,52,53] and policy trees [29]. Point-based techniques [38,42] approximate a POMDP's value function using a finite subset of beliefs which is iteratively updated. Algorithms include PBVI [38], Perseus [48], SARSOP [30] and HSVI [45]. Point-based methods can treat large POMDPs for both finiteand discounted infinite-horizon objectives [42].

Problem statement. In this paper we consider the problem: is the maximal expected total reward to reach a given goal state in a POMDP below a given threshold? We thus consider an infinite-horizon objective without discounting—also called an indefinite-horizon objective. A specific instance of the considered problem is the reachability probability to eventually reach a given goal state in a POMDP. This problem is undecidable [33,34] in general. Intuitively, this is due to the fact that POMDP policies need to consider the entire (infinite) observation history to make optimal decisions. For a POMDP, this notion is captured by an infinite, fully observable MDP, its belief MDP. This MDP is obtained from observation sequences inducing probabilities of being in certain states of the POMDP.

Previously proposed methods to solve the problem are e.g. to use approximate value iteration [22], optimisation and search techniques [1,12], dynamic programming [6], Monte Carlo simulation [43], game-based abstraction [51], and machine learning [13,14,19]. Other approaches restrict the memory size of the policies [35]. The synthesis of (possibly randomised) finite-memory policies is ETR-complete<sup>1</sup> [28]. Techniques to obtain finite-memory policies use e.g. parameter synthesis [28] or satisfiability checking and SMT solving [15,50].

Our approach. We tackle the aforementioned problem by computing underapproximations on maximal total expected rewards. This is done by considering finite unfoldings of the infinite belief MDP of the POMDP, and then applying abstraction. The key issue here is to find a suitable under-approximation of the POMDP's value function. We provide two techniques: a simple (cut-off) technique that uses a good policy on the POMDP, and a more advanced technique (belief clipping) that uses minimal shifts of probabilities between beliefs and can be applied on top of the simple approach. We use mixed-integer linear programming (MILP) to find such minimal probability shifts. Cut-off techniques for indefinite-horizon objectives have been used on computation trees—rather than on the belief MDP as used here—in Goal-HSVI [24]. Belief clipping amends the probabilities in a belief to be in a state of the POMDP yielding discretised values, i.e. an abstraction of the probability range [0, 1] is applied. Such grid-based approximations are inspired by Lovejoy's grid-based belief MDP discretisation method [32]. They have also been used in [7] in the context of dynamic programming for POMDPs, and to over-approximate the value function in model checking of POMDPs [8]. In fact, this paper on determining lower bounds for

<sup>1</sup> A decision problem is ETR-complete if it can be reduced to a polynomial-length sentence in the Existential Theory of the Reals (for which the satisfiability problem is decidable) in polynomial time, and there is such a reduction in the reverse direction.

indefinite-horizon objectives can be seen as the dual counterpart of [8]. Our key challenge—compared to the approach of [8]—is that the value at a certain belief cannot easily be under-approximated with a convex combination of values of nearby beliefs. On the other hand, an under-approximation can benefit from a "good" guess of some initial POMDP policy. In the context of [8], such a guessed policy is of limited use for over-approximating values in the POMDP induced by an optimal policy. Although our approach is applicable to all thresholds, the focus of our work is on determining under-approximations for quantitative objectives. Dedicated verification techniques for the qualitative setting—almost-sure reachability—are presented in [17,16,27].

Experimental results. We have implemented our cut-off and belief clipping approaches on top of the probabilistic model checker Storm [23] and applied it to a range of various benchmarks. We provide a comparison with the model checking approach in [37], and determine the tightness of our under-approximations by comparing them to over-approximations obtained using the algorithm from [8]. Our main findings from the experimental validation are:

# 2 Preliminaries and Problem Statement

Let Dist(A) := µ : A → [0, 1] | P <sup>a</sup>∈<sup>A</sup> µ(a) = 1 denote the set of probability distributions over a finite set A. The set supp(µ) := {a ∈ A | µ(a) > 0} is the support of µ ∈ Dist(A). Let R<sup>∞</sup> := R ∪ {∞, −∞}. We use Iverson bracket notation, where [x] = 1 if the Boolean expression x is true and [x] = 0 otherwise.

#### 2.1 Partially Observable MDPs

Definition 1 (MDP). A Markov decision process (MDP) is a tuple M = hS, Act, P, siniti with a (finite or infinite) set of states S, a finite set of actions Act, a transition function P: S ×Act ×S → [0, 1] with P s <sup>0</sup>∈<sup>S</sup> P(s, α, s<sup>0</sup> ) ∈ {0, 1} for all s ∈ S and α ∈ Act, and an initial state sinit.

We fix an MDP M := hS, Act, P, siniti. For s ∈ S and α ∈ Act, let post<sup>M</sup>(s, α) := {s <sup>0</sup> ∈ S | P(s, α, s<sup>0</sup> ) > 0} denote the set of α-successors of s in M. The set of enabled actions in s ∈ S is given by Act(s) := {α ∈ Act | post<sup>M</sup>(s, α) 6= ∅}.

Definition 2 (POMDP). A partially observable MDP (POMDP) is a tuple M = hM, Z, Oi, where M is the underlying MDP with |S| ∈ N, i.e. S is finite, Z is a finite set of observations, and O: S → Z is an observation function such that O(s) = O(s 0 ) =⇒ Act(s) = Act(s 0 ) for all s, s<sup>0</sup> ∈ S.

We fix a POMDP M := hM, Z, Oi with underlying MDP M. We lift the notion of enabled actions to observations z ∈ Z by setting Act(z) := Act(s) for some s ∈ S with O(s) = z which is valid since states with the same observations are required to have the same enabled actions. The notions defined for MDPs below also straightforwardly apply to POMDPs.

Remark 1. More general observation functions of the form O : S×Act → Dist(Z) can be encoded in this formalism by using a polynomially larger state space [16].

An infinite path through an MDP (and a POMDP) is a sequence π˜ = s0α1s1α<sup>2</sup> . . . such that αi+1 ∈ Act(si) and si+1 ∈ postM(s<sup>i</sup> , αi+1) for all i ∈ N. A finite path is a finite prefix πˆ = s0α<sup>1</sup> . . . αns<sup>n</sup> of an infinite path π˜. For finite πˆ let last(πˆ) := s<sup>n</sup> and |πˆ| := n. For infinite π˜ set |π˜| := ∞ and let π˜[i] denote the finite prefix of length i ∈ N. We denote the set of finite and infinite paths in M by Paths<sup>M</sup> fin and Paths<sup>M</sup> inf, respectively. Let Paths<sup>M</sup> := Paths<sup>M</sup> fin ∪ Paths<sup>M</sup> inf. Paths are lifted to the observation level by observation traces. The observation trace of a (finite or infinite) path π = s0α1s1α<sup>2</sup> . . . ∈ Paths<sup>M</sup> is O(π) := O(s0)α1O(s1)α<sup>2</sup> . . .. Two paths π, π<sup>0</sup> ∈ Paths<sup>M</sup> are observation-equivalent if O(π) = O(π 0 ).

Policies resolve the non-determinism present in MDPs (and POMDPs). Given a finite path πˆ, a policy determines the action to take at last(ˆπ).

Definition 3 (Policy). A policy for M is a function σ : Paths<sup>M</sup> fin → Dist(Act) such that for each path πˆ ∈ Paths<sup>M</sup> fin, supp(σ(ˆπ)) ⊆ Act(last(ˆπ)).

A policy σ is deterministic if |supp(σ(πˆ))| = 1 for all πˆ ∈ Paths<sup>M</sup> fin. Otherwise it is randomised. σ is memoryless if for all π, ˆ πˆ <sup>0</sup> ∈ Paths<sup>M</sup> fin we have last(πˆ) = last(πˆ 0 ) =⇒ σ(πˆ) = σ(πˆ 0 ). σ is observation-based if for all π, ˆ πˆ <sup>0</sup> ∈ Paths<sup>M</sup> fin it holds that O(πˆ) = O(πˆ 0 ) =⇒ σ(πˆ) = σ(πˆ 0 ). We denote the set of policies for M by Σ<sup>M</sup> and the set of observation-based policies for M by Σ<sup>M</sup> obs. A finite-memory policy (fm-policy) can be represented by a finite automaton where the current memory state and the state of the MDP determine the actions to take [4].

The probability measure µ σ,s <sup>M</sup> for paths in M under policy σ and initial state s is the probability measure of the Markov chain induced by M, σ, and s [4].

We use reward structures to model quantities like time, or energy consumption.

Definition 4 (Reward Structure). A reward structure for M is a function R: S × Act × S → R such that either for all s, s<sup>0</sup> ∈ S, α ∈ Act, R(s, α, s<sup>0</sup> ) ≥ 0 or for all s, s<sup>0</sup> ∈ S, α ∈ Act, R(s, α, s<sup>0</sup> ) ≤ 0 holds. In the former case, we call R positive, otherwise negative.

We fix a reward structure R for M. The total reward along a path π is defined as rewM,R(π) := P<sup>|</sup>π<sup>|</sup> <sup>i</sup>=1 R(si−1, α<sup>i</sup> , si). The total reward is always well-defined even if π is infinite—since all rewards are assumed to be either non-negative or non-positive. For an infinite path π˜ we define the total reward until reaching a set of goal states G ⊆ S by

$$\mathsf{rew}\_{M,\mathbf{R},G}(\tilde{\pi}) := \begin{cases} \mathsf{rew}\_{M,\mathbf{R}}(\hat{\pi}) & \text{if } \exists i \in \mathbb{N} : \hat{\pi} = \tilde{\pi}[i] \wedge last(\hat{\pi}) \in G \wedge \\ & \forall j < i : last(\tilde{\pi}[j]) \notin G, \\ \mathsf{rew}\_{M,\mathbf{R}}(\tilde{\pi}) & \text{otherwise}. \end{cases}$$

Intuitively, rewM,R,G(π˜) accumulates reward along π˜ until the first visit of a goal state s ∈ G. If no goal state is reached, reward is accumulated along the infinite path. The expected total reward until reaching G for policy σ and state s is

$$\mathsf{ER}\_{M,\mathbf{R}}^{\sigma}(s \mid = \Diamond G) := \int\_{\tilde{\pi} \in Paths\_{\mathrm{inf}}^{M}} \int\_{\mathrm{inf}} \mathsf{rew}\_{M,\mathbf{R},G}(\tilde{\pi}) \cdot \mu\_{M}^{\sigma,s}(d\tilde{\pi}).$$

Observation-based policies capture the notion that a decision procedure for a POMDP only accesses the observations and their history and not the entire state of the system. We are interested in reasoning about minimal and maximal values over all observation-based policies. For our explanations we focus on maximising (non-negative or non-positive) expected rewards. Minimisation can be achieved by negating all rewards.

Definition 5 (Maximal Expected Total Reward). The maximal expected total reward until reaching G from s in POMDP M is

$$\mathsf{ER}\_{\mathcal{M},\mathbf{R}}^{\max}(s \mathbb{H} \vDash \lozenge G) := \sup\_{\sigma \in \Sigma\_{\textup{obs}}^{\mathcal{M}}} \mathsf{ER}\_{\mathcal{M},\mathbf{R}}^{\sigma}(s \mathbb{H} \vDash \lozenge G).$$

We define ERmax <sup>M</sup>,R(♦G) := ERmax <sup>M</sup>,R(sinit |= ♦G).

The central problem of our work, the indefinite-horizon total reward problem, asks the question whether the maximal expected total reward until reaching a goal exceeds a given threshold.

Problem 1. Given a POMDP M, reward structure R, set of goal states G ⊆ S, and threshold λ ∈ R, decide whether ERmax <sup>M</sup>,R(♦G) ≤ λ.

Example 1. Fig. 1 shows a POMDP M with three states and two observations: O(s0) = O(s1) = and O(s2) = . A reward of 1 is collected when transitioning from s<sup>1</sup> to s<sup>2</sup> via the β-action. All other rewards are zero.

The policy that always selects α at s<sup>0</sup> and β at s<sup>1</sup> maximizes the expected total reward to reach G = {s2} but is not observation-based. The observation-based policy that for the first n ∈ N transition steps selects α and then selects β afterwards yields an expected total reward of 1 − ( <sup>1</sup>/2) <sup>n</sup>. With n → ∞ we obtain ERmax <sup>M</sup>,R(♦{s2}) = 1.

As computing maximal expected rewards exactly in POMDPs is undecidable [34], we aim at under-approximating the actual value ERmax <sup>M</sup>,R(♦G). This allows us to answer our problem negatively if the computed lower bound exceeds λ.

Remark 2.  Our approach can thus be used to obtain lower bounds on reachability probabilities in POMDPs. This also holds for almost-sure reachability (i.e. "is the reachability probabilty one?"), though dedicated methods like those presented in [17,16,27] are better suited for that setting.

#### 2.2 Beliefs

The semantics of a POMDP M are captured by its (fully observable) belief MDP. The infinite state space of this MDP consists of beliefs [3,44]. A belief is a distribution over the states of the POMDP where each component describes the likelihood to be in a POMDP state given a history of observations. We denote the set of all beliefs for M by B<sup>M</sup> := {b ∈ Dist(S) | ∀s, s<sup>0</sup> ∈ supp(b) : O(s) = O(s 0 )} and write O(b) ∈ Z for the unique observation O(s) of all s ∈ supp(b).

The belief MDP of M is constructed by starting in the belief corresponding to the initial state and computing successor beliefs to unfold the MDP. Let P(s, α, z) := P s <sup>0</sup>∈S [O(s 0 ) = z] · P(s, α, s<sup>0</sup> ) be the probability to observe z ∈ Z after taking action α in POMDP state s. Then, the probability to observe z after taking action α in belief b is P(b, α, z) := P s∈S b(s) · P(s, α, z). We refer to <sup>J</sup>b|α, z<sup>K</sup> ∈ BM—the belief after taking <sup>α</sup> in <sup>b</sup>, conditioned on observing <sup>z</sup>—as the α-z-successor of b. If P(b, α, z) > 0, it is defined component-wise as

$$[b|\alpha, z](s) := \frac{[O(s) = z] \cdot \sum\_{s' \in S} b(s') \cdot \mathbf{P}(s', \alpha, s)}{\mathbf{P}(b, \alpha, z)}$$

for all <sup>s</sup> <sup>∈</sup> <sup>S</sup>. Otherwise <sup>J</sup>b|α, z<sup>K</sup> is undefined.

Definition 6 (Belief MDP). The belief MDP of M is the MDP bel(M) = BM, Act, PB, binit , where B<sup>M</sup> is the set of all beliefs in M, Act is as for M, binit := {sinit 7→ 1} is the initial belief, and P<sup>B</sup> : B<sup>M</sup> × Act × B<sup>M</sup> → [0, 1] is the belief transition function with

$$\mathbf{P}^B(b,\alpha,b') := \begin{cases} \mathbf{P}(b,\alpha,z) & \text{if } b' = [b|\alpha, z],\\ 0 & \text{otherwise.} \end{cases}$$

We lift a POMDP reward structure R to the belief MDP [25].

Definition 7 (Belief Reward Structure). For beliefs b, b<sup>0</sup> ∈ B<sup>M</sup> and action α ∈ Act, the belief reward structure R<sup>B</sup> based on R associated with bel(M) is given by

$$\mathbf{R}^B(b,\alpha,b') := \frac{\sum\_{s \in S} b(s) \cdot \sum\_{s' \in S} [O(s') = O(b')] \cdot \mathbf{R}(s,\alpha,s') \cdot \mathbf{P}(s,\alpha,s')}{\mathbf{P}(b,\alpha,O(b'))}.$$

Given a set of goal states G ⊆ S, we assume—for simplicity—that there is a set of observations Z <sup>0</sup> ⊆ Z such that s ∈ G iff O(s) ∈ Z 0 . This assumption can always be ensured by transforming the POMDP M. See the full technical report [10] for details. The set of goal beliefs for G is given by G<sup>B</sup> := {b ∈ B<sup>M</sup> | supp(b) ⊆ G}.

We now lift the computation of expected rewards to the belief level. Based on the well-known Bellman equations [5], the belief MDP induces a function that maps every belief to the expected total reward accumulated from that belief.

Definition 8 (POMDP Value Function). For b ∈ BM, the n-step value function V<sup>n</sup> : B<sup>M</sup> → R of M is defined recursively as V0(b) := 0 and

$$V\_n(b) := [b \notin G\_{\mathcal{B}}] \cdot \max\_{\alpha \in Act} \sum\_{b' \in post^{\text{hol}(\mathcal{M})}(b,\alpha)} \mathbf{P}^B(b,\alpha,b') \cdot \left(\mathbf{R}^B(b,\alpha,b') + V\_{n-1}(b')\right).$$

Figure 2. Belief MDP bel(M) of POMDP M from Fig. 1

The (optimal) value function V ∗ : B<sup>M</sup> → R<sup>∞</sup> is given by V ∗ (b) := limn→∞ Vn(b).

The n-step value function is piecewise linear and convex [44]. Thus, the optimal value function can be approximated arbitrarily close by a piecewise linear convex function [47]. The value function yields expected total rewards in M and bel(M):

ERmax <sup>M</sup>,R(s |= ♦G) = ERmax bel(M),R<sup>B</sup> ({s 7→ 1} |= ♦GB) = V ∗ ({s 7→ 1}).

Example 2. Fig. 2 shows a fragment of the belief MDP of the POMDP from Fig. 1. Observe ERmax bel(M),R<sup>B</sup> (♦ {s<sup>2</sup> 7→ 1}) = 1.

We reformulate our problem statement to focus on the belief MDP.

Problem 2 (equivalent to Problem 1). For a POMDP M, reward structure R, goal states G ⊆ S, and threshold λ ∈ R, decide whether V ∗ ({sinit 7→ 1}) ≤ λ.

As the belief MDP is fully observable, standard results for MDPs apply. However, an exhaustive analysis of bel(M) is intractable since the belief MDP is—in general—infinitely large<sup>2</sup> .

### 3 Finite Exploration Under-Approximation

Instead of approximating values directly on the POMDP, we consider approximations of the corresponding belief MDP. The basic idea is to construct a finite abstraction of the belief MDP by unfolding parts of it and approximate values at beliefs where we decide not to explore. In the resulting finite MDP, under-approximative expected reward values can be computed by standard model checking techniques. We present two approaches for abstraction: belief cut-offs and belief clipping. We incorporate those techniques into an algorithmic framework that yields arbitrarily tight under-approximations.

The technical report [10] contains formal proofs of our claims.

<sup>2</sup> The set of all beliefs—i.e. the state space of bel(M)—is uncountable. The reachable fragment is countable, though, since each belief has at most |Z| many successors.

Figure 3. Applying belief cut-offs to the belief MDP from Fig. 2

#### 3.1 Belief Cut-Offs

The general idea of belief cut-offs is to stop exploring the belief MDP at certain beliefs—the cut-off beliefs—and assume that a goal state is immediately reached while sub-optimal reward is collected. Similar techniques have been discussed in the context of fully observable MDPs and other model types [11,26,49,2]. Our work adapts the idea of cut-offs for POMDP over-approximations described in [8] to under-approximations. The main idea of belief cut-offs shares similarities with the SARSOP [30] and Goal-HSVI [24] approaches. While they apply cut-offs on the level of the computation tree, our approach directly manipulates the belief MDP to yield a finite model.

Let V : B<sup>M</sup> → R<sup>∞</sup> with V(b) ≤ V ∗ (b) for all b ∈ BM. We call V an underapproximative value function and V(b) the cut-off value of b. In each of the cut-off beliefs b, instead of adding the regular transitions to its successors, we add a transition with probability 1 to a dedicated goal state bcut. In the modified reward structure R<sup>0</sup> , this cut-off transition is assigned a reward<sup>3</sup> of V(b), causing the value for a cut-off belief b in the modified MDP to coincide with V(b). Hence, the exact value of the cut-off belief—and thus the value of all other explored beliefs—is under-approximated.

Example 3. Fig. 3 shows the resulting finite MDP obtained when considering the belief MDP from Fig. 2 with single cut-off belief b = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ <sup>3</sup>/4}.

Computing cut-off values. The question of finding a suitable under-approximative value function V is central to the cut-off approach. For an effective approximation, such a function should be easy to compute while still providing values close to the optimum. If we assume a positive reward structure, the constant value 0 is always a valid under-approximation. A more sophisticated approach is to compute suboptimal expected reward values for the states of the POMDP using some arbitrary, fixed observation-based policy σ ∈ Σ<sup>M</sup> obs.  Then, we define the function U σ : B<sup>M</sup> → R<sup>∞</sup> as U σ (b) := P s∈supp(b) b(s) · U σ (s).

<sup>3</sup> We slightly deviate from Def. 4 by allowing transition rewards to be −∞ or +∞. Alternatively, we could introduce new sink states with a non-zero self-loop reward.

Lemma 1. U σ is an under-approximative value function, i.e. for all b ∈ BM:

$$\mathfrak{U}^{\sigma}(b) := \sum\_{s \in \operatorname{supp}(b)} b(s) \cdot U^{\sigma}(s) \le V^\*(b).$$

Thus, finding a suitable under-approximative value function reduces to finding "good" policies for M, e.g. by using randomly guessed fm-policies, machine learning methods [13], or a transformation to a parametric model [28].

#### 3.2 Belief Clipping

The cut-off approach provides a universal way to construct an MDP which underapproximates the expected total reward value for a given POMDP. The quality of the approximation, however, is highly dependent on the under-approximative value function used. Furthermore, regions where the belief MDP slowly converges towards a belief may pose problems in practice.

As a potential remedy for these problems, we propose a different concept called belief clipping. Intuitively, the procedure shifts some of the probability mass of a belief b in order to transform b to another belief ˜b. We then connect b to ˜b in a way that the accuracy of our approximation of the value V ∗ (b) depends only on the approximation of V ∗ ( ˜b) and the so-called clipping value—some notion of distance between b and ˜b that we discuss below. We can thus focus on exploring the successors of ˜b to obtain good approximations for both beliefs b and ˜b.

Definition 9 (Belief Clip). For b ∈ BM, we call µ: supp(b) → [0, 1] a belief clip if ∀s ∈ supp(b): µ(s) ≤ b(s) and P(µ) := P <sup>s</sup>∈supp(b) µ(s) < 1. The belief (b µ) ∈ B<sup>M</sup> induced by µ is defined by

$$\forall s \in \operatorname{supp}(b) \colon \ (b \ominus \mu)(s) \ := \ \frac{b(s) - \mu(s)}{1 - \sum(\mu)}.$$

Intuitively, a belief clip µ for b describes for each s ∈ supp(b) the probability mass that is removed ("clipped away") from b(s). The induced belief is obtained when normalising the resulting values so that they sum up to one.

Example 4. For belief b = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ <sup>3</sup>/4}, consider the two belief clips µ<sup>1</sup> = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ <sup>1</sup>/4} and µ<sup>2</sup> = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ 0}. Both induce the same belief: (b µ1) = (b µ2) = {s<sup>0</sup> 7→ 0, s<sup>1</sup> 7→ 1}.

We have supp((b µ)) ⊆ supp(b), which also implies O((b µ)) = O(b). Given some candidate belief ˜b, consider the set of inducing belief clips:

$$\mathcal{L}(b,\tilde{b}) := \left\{ \mu \colon \operatorname{supp}(b) \to [0,1] \mid \mu \text{ is a belief clip for } b \text{ with } \tilde{b} = (b \ominus \mu) \right\}.$$

Belief ˜b is called an adequate clipping candidate for b iff C(b, ˜b) 6= ∅.

Definition 10 (Clipping Value). For b ∈ B<sup>M</sup> and adequate clipping candidate ˜b, the clipping value is <sup>∆</sup>b→˜<sup>b</sup> := P(δb→˜<sup>b</sup> ), where δb→˜<sup>b</sup> := arg minµ∈C(b,˜b) P(µ). The values δb→˜<sup>b</sup> (s) for s ∈ supp(b) are the state clipping values.

Figure 4. Applying belief clipping to the belief MDP from Fig. 2

Given a belief b and an adequate clipping candidate ˜b, we outline how the notion of belief clipping is used to obtain valid under-approximations. We assume b 6= ˜b, implying 0 < ∆b→˜<sup>b</sup> < 1. Instead of exploring all successors of b in bel(M), the approach is to add a transition from b to ˜b. The newly added transition has probability 1 − ∆b→˜<sup>b</sup> and gets assigned a reward of 0. The remaining probability mass (i.e. ∆b→˜<sup>b</sup> ) leads to a designated goal state bcut. To guarantee that—in general—the clipping procedure yields a valid under-approximation, we need to add a corrective reward value to the transition from b to bcut. Let L : S → R<sup>∞</sup> which maps each POMDP state to its minimum expected reward in the underlying, fully observable MDP M of M<sup>4</sup> , i.e. L(s) = ERmin M,R(s |= ♦G). This function soundly under-approximates the state values which can be achieved by any observation-based policy. It can be generated using standard MDP analysis. Given state clipping values δb→˜<sup>b</sup> (s) for s ∈ supp(b), the reward for the transition from b to bcut is P s∈supp(b) (δb→˜<sup>b</sup> (s)/∆b→˜<sup>b</sup> ) · L(s).

Example 5. For the belief MDP from Fig. 2, belief b = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ <sup>3</sup>/4}, and clipping candidate ˜<sup>b</sup> <sup>=</sup> {s<sup>0</sup> 7→ <sup>0</sup>, s<sup>1</sup> 7→ <sup>1</sup>} we get <sup>∆</sup>b→˜<sup>b</sup> <sup>=</sup> <sup>1</sup>/4, as <sup>δ</sup>b→˜<sup>b</sup> <sup>=</sup> µ<sup>2</sup> = {s<sup>0</sup> 7→ <sup>1</sup>/4, s<sup>1</sup> 7→ 0} with the belief clip µ<sup>2</sup> as in Example 4. Furthermore, L(s0) = 0. The resulting MDP following our construction above is given in Fig. 4.

The following lemma shows that the construction yields an under-approximation.

$$\mathbf{Lemma 2. } (1 - \Delta\_{b \to \tilde{b}}) \cdot V^\*(\tilde{b}) \, + \, \Delta\_{b \to \tilde{b}} \cdot \sum\_{s \in supp(b)} \frac{\delta\_{b \to \tilde{b}}(s)}{\Delta\_{b \to \tilde{b}}} \cdot \mathfrak{L}(s) \, \le \, \, V^\*(b).$$

Proof (sketch). To gain some intuition, consider the special case, where ∆b→˜<sup>b</sup> = δb→˜<sup>b</sup> (s) = b(s) for some s ∈ supp(b). The clipping candidate ˜b can be interpreted as the conditional probability distribution arising from distribution b given that s is not the current state. The value V ∗ (b) can be split into the sum of (i) the probability that s is not the current state times the reward accumulated from belief ˜b and (ii) the probability that s is the current state times the reward accumulated from s, i.e. from the belief {s 7→ 1}. However, for the two summands

<sup>4</sup> When rewards are negative, we might have L(s) = −∞ for many s ∈ S \ G in which case the applicability of the clipping approach is very limited.

∀b

we must consider a policy that does not distinguish between the beliefs b, ˜b, and {s 7→ 1} as well as their observation-equivalent successors. In other words, the same sequence of actions must be executed when the same observations are made.

We consider such a policy that in addition is optimal at ˜b, i.e. the reward accumulated from ˜b is equal to V ∗ ( ˜b). For the reward accumulated from {s 7→ 1}, L(s) provides a lower bound. Hence, (1 − b(s)) · V ∗ (b) + b(s) · L(s) is a lower bound for the reward accumulated from b. A formal proof is given in [10]. ut

To find a suitable clipping candidate for a given belief b, we consider a finite candidate set B ⊆ B<sup>M</sup> consisting of beliefs with observation O(b). These beliefs do not need to be reachable in the belief MDP. The set can be constructed, e.g. by taking already explored beliefs or by using a fixed, discretised set of beliefs.

We are interested in minimising the clipping value ∆b→<sup>b</sup> <sup>0</sup> over all candidate beliefs b <sup>0</sup> ∈ B. A naive approach is to explicitly compute all clipping values for all candidates. We are using mixed-integer linear programming (MILP) [41] instead. An MILP is a system of linear inequalities (constraints) and a linear objective function considering real-valued and integer variables. A feasible solution of the MILP is a variable assignment that satisfies all constraints. An optimal solution is a feasible solution that minimises the objective function.

Definition 11 (Belief Clipping MILP). The belief clipping MILP for belief b ∈ B<sup>M</sup> and finite set of candidates B ⊆ {b <sup>0</sup> ∈ B<sup>M</sup> | O(b 0 ) = O(b)} is given by:

minimise ∆ such that: X b <sup>0</sup>∈B ab <sup>0</sup> = 1 . Select exactly one candidate b 0 (1)

<sup>0</sup> ∈ B: a<sup>b</sup> <sup>X</sup> <sup>0</sup> ∈ {0, <sup>1</sup>} (2)

$$\sum\_{s \in supp(b)} \delta\_s = \Delta \qquad \Rightarrow Compute \ clipping \ value \ for \ selected \ b' \tag{3}$$

$$\begin{array}{ccc} \forall s \in support(b) \colon & \delta\_s \in [0, b(s)] \\ \forall b' \in \mathfrak{B} \colon & \delta\_s \ge b(s) - (1 - \Delta) \cdot b'(s) - (1 - a\_{b'}) \end{array} \tag{4}$$

The MILP consists of O(|supp(b)| + |B|) variables and O(|supp(b)| · |B|) constraints. For b <sup>0</sup> ∈ B, the binary variable a<sup>b</sup> <sup>0</sup> indicates whether b <sup>0</sup> has been chosen as the clipping candidate. Moreover, we have variables δ<sup>s</sup> for s ∈ supp(b) and a variable ∆ to represent the (state) clipping values for b and the chosen candidate b 0 . Constraints 1 and 2 enforce that exactly one of the a<sup>b</sup> <sup>0</sup> variables is one, i.e. exactly one belief is chosen. Constraint 3 forces ∆ to be the sum of all state clipping values. δ<sup>s</sup> variables get a value between zero and b(s) (Constraint 4). Constraint 5 only affects δ<sup>s</sup> if the corresponding belief is chosen. Otherwise, a<sup>b</sup> 0 is set to 0 and the value on the right-hand side becomes negative. If a belief b 0 is chosen, the minimisation forces Constraint 5 to hold with equality as the right-hand side is greater or equal to 0. Assuming ∆ is set to a value below 1, we obtain a valid clipping values as

$$\forall s \in support(b) \colon \quad \delta\_s = b(s) - (1 - \Delta) \cdot b'(s) \quad \Longleftrightarrow \quad b'(s) = \frac{b(s) - \delta\_s}{1 - \Delta}.$$

Input : POMDP M = hM, Z, Oi with M = hS, Act, P, siniti, reward structure R, goal states G ⊆ S, under-approx. value function V, function L : S → R <sup>∞</sup> with L(s) = ERmin M,R(s |= ♦G) Output : Clipping belief MDP K<sup>M</sup> and reward structure R<sup>K</sup> 1 S <sup>K</sup> ← {binit, bcut} with binit = {sinit 7→ 1} and a new belief state bcut 2 P <sup>K</sup>(bcut, cut, bcut) ← 1, R<sup>K</sup>(bcut, cut, bcut) ← 0 // add self-loop 3 Q ← {binit} // initialize exploration set 4 while Q 6= ∅ do 5 b ← chooseBelief(Q), Q ← Q \ {b} // pop next belief to explore from Q 6 if supp(b) ⊆ G then P <sup>K</sup>(b, goal, b) ← 1, R<sup>K</sup>(b, goal, b) ← 0 // add self-loop 7 else if exploreBelief(b) then // expand b <sup>8</sup> foreach α ∈ Act(b) do // Using bel(M) and R<sup>B</sup> as in Defs. 6 and 7 9 foreach b <sup>0</sup> ∈ post bel(M) (b, α) do 10 P <sup>K</sup>(b, α, b<sup>0</sup> ) ← P <sup>B</sup>(b, α, b<sup>0</sup> ), R<sup>K</sup>(b, α, b<sup>0</sup> ) ← R<sup>B</sup>(b, α, b<sup>0</sup> ) 11 if b <sup>0</sup> ∈/ S <sup>K</sup> then S <sup>K</sup> ← S <sup>K</sup> ∪ {b 0 }, Q ← Q ∪ {b 0 } 12 else // apply cut-off and clipping to b 13 P <sup>K</sup>(b, cut, bcut) ← 1, R<sup>K</sup>(b, cut, bcut) ← V(b) // add cut-off transition 14 choose a finite set B ⊆ B<sup>M</sup> of clipping candidates for b <sup>15</sup> ˜b, ∆b→˜<sup>b</sup> , δb→˜<sup>b</sup> ← solveClippingMILP(b, B) <sup>16</sup> if ˜b 6= b and ˜b is adequate then // Clip b using ˜b 17 P <sup>K</sup>(b, clip, ˜b) <sup>←</sup> (1−∆b→˜<sup>b</sup> ), P <sup>K</sup>(b, clip, bcut) ← ∆b→˜<sup>b</sup> <sup>18</sup> R<sup>K</sup>(b, clip, ˜b) ← 0, R<sup>K</sup>(b, clip, bcut) ← P s∈supp(b) δ b→˜b (s) ∆b→˜<sup>b</sup> · L(s) <sup>19</sup> if ˜b /∈ S <sup>K</sup> then S <sup>K</sup> ← S <sup>K</sup> ∪ {˜b}, Q ← Q ∪ {˜b}

20 return K<sup>M</sup> = S <sup>K</sup>, Act ] {goal, cut, clip} , P <sup>K</sup>, binit and R<sup>K</sup>

Algorithm 1: Belief exploration algorithm with cut-offs and clipping

A trivial solution of the MILP is always obtained by setting a<sup>b</sup> <sup>0</sup> and ∆ to 1 and δ<sup>s</sup> to b(s) for all s and an arbitrary b <sup>0</sup> ∈ B. This corresponds to an invalid belief clip. However, as we minimise the value for ∆, we can conclude that no belief in the candidate set is adequate for clipping if ∆ is 1 in an optimal solution.

Theorem 1. An optimal solution to the belief clipping MILP for belief b and candidate set B sets a˜<sup>b</sup> to 1 and ∆ to a value below 1 iff ˜b ∈ B is an adequate clipping candidate for b with minimal clipping value.

#### 3.3 Algorithm

We incorporate belief cut-offs and belief clipping into an algorithmic framework outlined in Algorithm 1. As input, the algorithm takes an instance of Problems 1 and 2, i.e. a POMDP M with reward structure R and goal states G. In addition, the algorithm considers an under-approximative value function V (Sect. 3.1) and a function L for the computation of corrective reward values (Sect. 3.2).

Lines 1 and 2 initialise the state set S <sup>K</sup> of the under-approximative MDP K<sup>M</sup> with the initial belief binit and the designated goal state bcut which has only one transition to itself with reward 0. Furthermore, we initialise the exploration set Q by adding binit (Line 3). During the computation, Q is used to keep track of all beliefs we still need to process. We then execute the exploration loop (Lines 4 to 19) until Q becomes empty. In each exploration step, a belief b is selected<sup>5</sup> and removed from Q. There are three cases for the currently processed belief b.

If supp(b) ⊆ G, i.e. b is a goal belief, we add a self-loop with reward 0 to b and continue with the next belief (Line 6). b is not expanded as successors of goal beliefs will not influence the result of the computation.

If b is not a goal belief, we use a heuristic function<sup>6</sup> exploreBelief to decide if b is expanded in Line 7. Lines 8 to 11 outline the expansion step. The transitions from b to its successor beliefs and the corresponding rewards as in the original belief MDP (see Sect. 2.2) are added. Furthermore, the successor beliefs that have not been encountered before are added to the set of states S <sup>K</sup> and the exploration set Q.

If b is not expanded, we apply the cut-off approach and the clipping approach to b in Lines 12 to 19. In Line 13 we add a cut-off transition from b to bcut with a new action cut. We use the given under-approximative value function V to compute the cut-off reward. Towards the clipping approach, a set of candidate beliefs is chosen and the belief clipping MILP for b and the candidate set is constructed as described in Def. 11 (Lines 14 and 15). If an adequate candidate ˜b with clipping values ∆b→˜<sup>b</sup> and δb→˜<sup>b</sup> (s) for s ∈ supp(b) has been found, we add the transitions from b to bcut and to ˜b using a new action clip and probabilities ∆b→˜<sup>b</sup> and 1 − ∆b→˜<sup>b</sup> , respectively. Furthermore, we equip the transitions with reward values as described in Sect. 3.2 using the given function L (Lines 16 to 18). If the clipping candidate ˜b has not been encountered before, we add it to the state space of the MDP and to the exploration set in Line 19.

The result of the algorithm is an MDP K<sup>M</sup> with reward structure R<sup>K</sup>. The set of states S <sup>K</sup> of K<sup>M</sup> contains all encountered beliefs. To guarantee termination of the algorithm, the decision heuristic exploreBelief has to stop exploring further beliefs at some point. Moreover, the handling of clipping candidates in Line 19 should not add new beliefs to Q infinitely often. We therefore fix a finite set of candidate beliefs B # ⊆ B<sup>M</sup> and make sure that the candidate sets B in Line 14 satisfy (B \ S <sup>K</sup>) ⊆ B#. To ensure a certain progress in the exploration "clip-cycles"—i.e. paths of the form b<sup>1</sup> clip . . . clip b<sup>n</sup> clip b1—are avoided in KM. This can be done, e.g. by always expanding the candidate beliefs b ∈ B#.

Expected total rewards until reaching the extended set of goal beliefs Gcut := G<sup>B</sup> ∪ {bcut} in K<sup>M</sup> under-approximate the values in the belief MDP:

Theorem 2. For all beliefs b ∈ S <sup>K</sup> \ {bcut} it holds that

$$\mathsf{ER}\_{\mathsf{K},\mathcal{M},\mathbf{R}^{\mathcal{K}}}^{\max}(b \models \Diamond G\_{cut}) \leq V^\*(b) = \mathsf{ER}\_{bel(\mathcal{M}),\mathbf{R}^B}^{\max}(b \models \Diamond G\_{\mathcal{B}}).$$

Corollary 1. ERmax <sup>K</sup>M,R<sup>K</sup> (♦Gcut) <sup>≤</sup> ERmax <sup>M</sup>,R(♦G).

<sup>5</sup> For example, Q can be implemented as a FIFO queue.

<sup>6</sup> The decision can be made for example by considering the size of the already explored state space such that the expansion is stopped if a size threshold has been reached. More involved decision heuristics are subject to further research.

Table 1. Results for benchmark POMDPs with maximisation objective

#### 4 Experimental Evaluation

Implementation details. We integrated Algorithm 1 in the probabilistic model checker Storm [23] as an extension of the POMDP verification framework described in [8]. Inputs are a POMDP—encoded either explicitly or using an extension of the Prism language [37]—and a property specification. Internally, POMDPs and MDPs are represented using sparse matrices. The implementation supports minimisation<sup>7</sup> and maximisation of reachability probabilities, reachavoid probabilities (i.e. the probability to avoid a set of bad state until a set of goal states is reached), and expected total rewards. In a preprocessing step, functions V and L as considered in Algorithm 1 are generated. For V, we consider the function U <sup>σ</sup> as in Lemma 1, where σ is a memoryless observation-based policy given by a heuristic<sup>8</sup> . For the function L, we apply standard MDP analysis on the underlying MDP. When exploring the abstraction MDP KM, our heuristic expands a belief iff |S <sup>K</sup>| ≤ |S|·maxz∈<sup>z</sup> |O<sup>−</sup><sup>1</sup> (z)|, where |S <sup>K</sup>| is the number of already explored beliefs and |O<sup>−</sup><sup>1</sup> (z)| is the number of POMDP states with observation z. Belief clipping can either be disabled entirely, or we consider candidate sets B ⊆ B# η , where B # η := {b ∈ B | ∀s ∈ S : b(s) ∈ {i/<sup>η</sup> | i ∈ N, 0 ≤ i ≤ η}} forms a finite, regular grid of beliefs with resolution η ∈ N \ {0}. Grid beliefs b ∈ B# <sup>η</sup> are always expanded.

<sup>7</sup> For minimisation, the under-approximation yields upper bounds.

<sup>8</sup> The heuristic uses optimal values obtained on the fully observable underlying MDP.

Table 2. Results for benchmark POMDPs with minimisation objective

Furthermore, we exclude clipping candidates ˜<sup>b</sup> with <sup>δ</sup>b→˜<sup>b</sup> (s) > 0 for s with L(s) = −∞; clipping with such candidates is not useful as it induces a value of −∞. Expected total rewards on fully observable MDPs are computed using Sound Value Iteration [39] with relative precision 10−<sup>6</sup> . MILPs are solved using Gurobi [21].

Set-up. We evaluate our under-approximation approach with cut-offs only and with enabled belief clipping procedure using grid resolutions η = 2, 3, 4, 6. We consider the same POMDP benchmarks<sup>9</sup> as in [37,8]. The POMDPs are scalable versions of case studies stemming from various application domains. To establish an external baseline, we compare with the approach of [37] implemented in Prism [31]. Prism generates an under-approximation based on an optimal policy for an over-approximative MDP which—in contrast to Storm—means that always both, under- and over-approximations, have to be computed. We ran Prism with resolutions η = 2, 3, 4, 6, 8, 10 and report on the best approximation obtained. To provide a further reference for the tightness of our under-approximation, we compute over-approximative bounds as in [8] using the implementation in Storm with a resolution of η = 8. All experiments were run on an Intel® Xeon® Platinum 8160 CPU using 4 threads<sup>10</sup>, 64GB RAM and a time limit of 2 hours.

Results. Tables 1 and 2 show our results for maximising and minimising properties, respectively. The first columns contain for each POMDP the benchmark name,

<sup>9</sup> Instances with a finite belief MDP that would be fully explored by our algorithm are omitted since the exact value can be obtained without approximation techniques.

<sup>10</sup> For our implementation, only Gurobi runs multi-threaded. Prism uses multiple threads for garbage collection.

Figure 5. Accuracy for Drone 4-2 with different sizes of approximation MDP K<sup>M</sup>

model parameters, property type (probabilities (P) or rewards (R)), and the numbers of states, state-action pairs, and observations. Column Prism gives the result with the smallest gap between over- and under-approximation computed with the approach of [37]. For maximising (minimising) properties, our approach competes with the lower (upper) bound of the provided interval. The relevant value is marked in bold. We also provide the computation time and the considered resolution η. For our implementation, we give results for the configuration with disabled clipping and for clipping with different resolutions η. In each cell, we give the obtained value, the computation time and the number of states in the abstraction MDP KM. Time- and memory-outs are indicated by TO and MO. The right-most column indicates the over-approximation value computed via [8]. Discussion. The pure cut-off approach yields valid under-approximations in all benchmark instances—often exceeding the accuracy of the approach of [37] while being consistently faster. In some cases, the resulting values improve when clipping is enabled. However, larger candidate sets significantly increase the computation time which stems from the fact that many clipping MILPs have to be solved.

For Drone 4-2, Fig. 5 plots the resulting under-approximation values (y-axis) for varying sizes of the explored MDP K<sup>M</sup> (x-axis). The horizontal, dashed line indicates the computed over-approximation value. The quality of the approximation further improves with an increased number of explored beliefs.

#### 5 Conclusion

We presented techniques to safely under-approximate expected total rewards in POMDPs. The approach scales to large POMDPs and often produces tight lower bounds. Belief clipping generally does not improve on the simpler cut-off approach in terms of results and performance. However, considering—and optimising—the approach for particular classes of POMDPs might prove beneficial. Future work includes integrating the algorithm into a refinement loop that also considers over-approximation techniques from [8]. Furthermore, lifting our approach to partially observable stochastic games is promising.

Data Availability. The artifact [9] accompanying this paper contains source code, benchmark files, and replication scripts for our experiments.

### References

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# Correct Probabilistic Model Checking with Floating-Point Arithmetic?

Arnd Hartmanns

University of Twente, Enschede, The Netherlands<EMAIL_ADDRESS>
Abstract. Probabilistic model checking computes probabilities and expected values related to designated behaviours of interest in Markov models. As a formal verification approach, it is applied to critical systems; thus we trust that probabilistic model checkers deliver correct results. To achieve scalability and performance, however, these tools use finite-precision floating-point numbers to represent and calculate probabilities and other values. As a consequence, their results are affected by rounding errors that may accumulate and interact in hard-to-predict ways. In this paper, we show how to implement fast and correct probabilistic model checking by exploiting the ability of current hardware to control the direction of rounding in floating-point calculations. We outline the complications in achieving correct rounding from higherlevel programming languages, describe our implementation as part of the Modest Toolset's mcsta model checker, and exemplify the tradeoffs between performance and correctness in an extensive experimental evaluation across different operating systems and CPU architectures.

# 1 Introduction

Given a Markov chain or Markov decision process (MDP [25]) model of a safetyor performance-critical system, probabilistic model checking (PMC) calculates quantitative properties of interest: the probability of (rare or catastrophic) failures, the expected recovery time after service interruption, or the long-run average throughput. These properties involve probabilities or expected costs/rewards of sets of model behaviours, and are often specified in a temporal logic like PCTL [16]. As a formal verification approach, users place great trust in the results delivered by a PMC tool such as Prism [22], Storm [9], ePMC [15], or the Modest Toolset's [18] mcsta. In contrast to classical model checkers for functional, Boolean-valued properties specified in e.g. LTL or CTL [2], a probabilistic model checker is inherently quantitative: the input model contains real-valued probabilities and costs/rewards; PCTL makes comparisons between real-valued constants and probabilities; the most efficient algorithms numerically iterate towards a fixpoint; and the final result itself may well be a real number.

<sup>?</sup> This work was supported by NWO VENI grant no. 639.021.754 and the EU's Horizon 2020 research and innovation programme under MSCA grant agreement 101008233.

Often, we can restrict to rationals, which simplifies the theory and facilitates "exact" algorithms using arbitrary-precision rational number datatypes. These algorithms only work for small models (as shown in the most recent QComp 2020 competition of quantitative verification tools [6]). In this paper, we thus focus on the PMC techniques that scale to large problems: those building upon iterative numerical algorithms, in particular value iteration (VI) [8]. We restrict to probabilistic reachability, i.e. calculating the probability to eventually reach a goal state, as this is the core problem in PMC for MDP. Embedded in the usual recursive CTL algorithm, it allows us to check any (unbounded) PCTL formula.

Starting from a trivial underapproximation of the reachability probability for each state of the model, VI iteratively improves the value of each state based on its successors' values. The true reachability probabilities are the least fixpoint of this procedure, towards which the algorithm converges. For roughly a decade, PMC tools implemented VI by stopping once the relative or absolute difference between subsequent iterations was below a threshold . Haddad and Monmege [12] showed in 2014<sup>1</sup> that this does not guarantee a difference of ≤ between the reported and the true probability, putting in question the trust placed in PMC tools. Then variants of VI were developed that provide sound, i.e. -correct, results: interval iteration (II) [3,5,13], sound value iteration (SVI) [26], and optimistic value iteration (OVI) [19]. We focus on II as the prototypical sound algorithm. It additionally iterates on an overapproximation; its stopping criterion is the difference between over- and underapproximation being ≤ .

If all probabilities in an MDP are rational numbers, then the true reachability probability as well as all intermediate values in II are rational, too. Yet implementing II with arbitrary-precision rationals is impractical since the smallerand-smaller differences between intermediate values end up using excessive computation time and memory. II is thus implemented with fixed-precision (usually 64-bit IEEE 754 double precision) floating point numbers. These, however, cannot represent all rationals, so operations must round to nearby representable values. Although II is numerically benign, consisting only of multiplications and additions within [0, 1], the default round to nearest, ties to even policy can cause II to deliver incorrect results. Wimmer et al. [29] show an example where PMC tools incorrectly state that a simple PCTL property is satisfied by a small Markov chain due to the underlying numeric difference having disappeared in rounding. We confirmed with current versions of Prism, Storm, and mcsta that the problem persists to today, even when requesting a "sound" algorithm like II. Wimmer et al. propose interval arithmetic to avoid such problems, cautioning that

[...] the memory consumption will roughly double, since two numbers for the interval bounds have to be stored [...]. The runtime will be higher by a small factor, because we need to derive lower and upper bounds for the intervals, requiring two model checking runs per sub-formula. [29, p. 5]

They did not provide an implementation, and we are not aware of any to date.

<sup>1</sup> Wimmer et al. [29] already in 2008 mention this problem in a more general setting, but neither give a concrete counterexample nor propose a solution tailored to PMC.

Our contribution. We present the first PMC implementation that computes correct lower and upper bounds on reachability probabilities despite using floatingpoint arithmetic. We benefit from two developments since Wimmer et al.'s paper of 2008: First, II (published 2014) already uses intervals (though not as Wimmer et al. envisioned), necessarily doubling memory consumption compared to VI (as do SVI and OVI, so it appears an unavoidable cost of soundness). In place of "two model checking runs per sub-formula", we can make the two interleaved computations inside II safe w.r.t. rounding. Second, hardware and programming language support for controlling the rounding direction in floating-point operations has improved, in particular with the AVX-512 instruction set in the newest x86-64 CPUs and widespread compiler support for C99's "floating-point environment" header fenv.h. Nevertheless, it is nontrivial to achieve runtime that is only "higher by a small factor". For the analysis of probabilistic systems, the only related use of safe rounding we are aware of is in the SSMT tool SiSAT [27].

Structure. We recap PMC and II (Sect. 2) as well as problems and solutions related to rounding in floating-point arithmetic in Sect. 3. We then present our new approach in Sect. 4, including important implementation aspects. The performance of our approach is crucial to its adoption in tools; thus in Sect. 5 we report on extensive experiments across different software and hardware configurations on models from the Quantitative Verification Benchmark Set (QVBS) [20].

#### 2 Probabilistic Model Checking

We write { x<sup>1</sup> 7→ y1, . . . } to denote the function that maps all x<sup>i</sup> to y<sup>i</sup> . Given a set S, its powerset is 2 <sup>S</sup>. A (discrete) probability distribution over S is a function µ ∈ S → [0, 1] with countable support spt(µ) def <sup>P</sup> <sup>=</sup> { <sup>s</sup> <sup>∈</sup> <sup>S</sup> <sup>|</sup> <sup>µ</sup>(s) <sup>&</sup>gt; <sup>0</sup> } and <sup>s</sup>∈spt(µ) µ(s) = 1. Dist(S) is the set of all probability distributions over S. If µ(s) ∈ Q for all s ∈ S, we call µ a rational probability distribution, in DistQ(S).

Markov decision processes (MDP) [25] combine the nondeterminism of Kripke structures with the finite random choices of discrete-time Markov chains (DTMC).

Definition 1. A Markov decision process (MDP) is a triple M = hS, s<sup>I</sup> , Ti where S is a finite set of states with initial state s<sup>I</sup> ∈ S and T : S → 2 DistQ(S) is the transition function. T(s) must be finite and non-empty for all s ∈ S.

For s ∈ S, an element µ of T(s) is a transition, and if s <sup>0</sup> ∈ spt(µ), then the transition has a branch to successor state s <sup>0</sup> with probability µ(s 0 ). If |T(s)| = 1 for all s ∈ S, then M is a DTMC.

Example 1. Fig. 1 shows our example MDP M<sup>γ</sup> n , which is actually a DTMC. It is a simplified and parametrised version of the counterexample of Wimmer et al. [29, Fig. 2]. It is parametrised in terms of n ∈ N (determining the number of chained states with transitions labelled b) and γ ∈ (0, 0.5) (changing some probabilities). We draw transitions as lines to an intermediate node from which

Fig. 1. Example parametrised MDP M<sup>γ</sup> n

probability-labelled branches lead to successor states. We omit the intermediate node for transitions with a single branch, and label some transitions to easily refer to them. M<sup>γ</sup> <sup>n</sup> has 4 + n states and transitions, and 7 + 2n branches.

In practice, higher-level modelling languages like Modest [14] are used to specify MDP. The semantics of an MDP is captured by its paths. A path represents a concrete resolution of all nondeterministic and probabilistic choices. Formally:

Definition 2. A finite path is a sequence πfin = s<sup>0</sup> µ<sup>0</sup> s<sup>1</sup> µ<sup>1</sup> . . . µn−1s<sup>n</sup> where s<sup>i</sup> ∈ S for all i ∈ { 0, . . . , n } and µ<sup>i</sup> ∈ T(si)∧µi(si+1) > 0 for all i ∈ { 0, . . . , n− 1 }. Let |πfin| def = n and last(πfin) def = sn. Πfin(s) is the set of all finite paths starting in s. A path is an analogous infinite sequence π, and Π(s) is the set of all paths starting in s. We write s ∈ π if ∃ i: s = si.

A scheduler (or adversary, policy or strategy) only resolves the nondeterministic choices of M. For this paper, memoryless deterministic schedulers suffice [4].

Definition 3. A function s : S → Dist(S) is a scheduler if, for all s ∈ S, we have s(s) ∈ T(s). The set of all schedulers of M is S(M).

We are interested in reachability probabilities. Let M|<sup>s</sup> = hS, s<sup>I</sup> , T|si with T|s(s) = { s(s) } be the DTMC induced by s on M. Via the standard cylinder set construction [10, Sect. 2.2] on M|s, a scheduler induces probability measures P M,s s on measurable sets of paths starting in s ∈ S.

Definition 4. For state s and goal state g ∈ S, the maximum and minimum probability of reaching g from s is defined as PM,s max ( g) = sups∈<sup>S</sup> P M,s <sup>s</sup> ({ π ∈ Π(s) | g ∈ π }) and P M,s min ( <sup>g</sup>) = infs∈<sup>S</sup> <sup>P</sup> M,s <sup>s</sup> ({ π ∈ Π(s) | g ∈ π }), respectively.

The definition extends to sets G of goal states. We omit the superscript for M when it is clear from the context, and if we omit that for s, then s = s<sup>I</sup> . From now on, whenever we have an MDP with a set of goal states G, we assume w.l.o.g. that all g ∈ G are absorbing, i.e. every g only has one self-loop transition.

Definition 5. A maximal end component (MEC) of M is a maximal (sub-)MDP hS 0 , T<sup>0</sup> , s<sup>0</sup> I i where S <sup>0</sup> ⊆ S, T 0 (s) ⊆ T(s) for all s ∈ S 0 , and the directed graph with vertex set S <sup>0</sup> and edge set { hs, s<sup>0</sup> i | ∃ µ ∈ T 0 (s): µ(s 0 ) > 0 } is strongly connected.

```
1 function II(M = hS, sI , Ti, G, opt, )
      // Preprocessing
 2 if opt = max then M := CollapseMECs(M, G) // collapse MECs
 3 S0 := Prob0(M, G, opt), S1 := Prob1(M, G, opt) // identify 0/1 states
 4 l := { s 7→ 0 | s ∈ S \ S1 } ∪ { s 7→ 1 | s ∈ S1 } // initialise lower vector
 5 u := { s 7→ 0 | s ∈ S0 } ∪ { s 7→ 1 | s ∈ S \ S0 } // initialise upper vector
      // Iteration
 6 while (u(sI ) − l(sI ))/l(sI ) >  do // while relative error > :
 7 foreach s ∈ S \ (S0 ∪ S1) do // update non-0/1 states:
 8 l(s) := optµ∈T (s)
                            P
                              s0∈spt(µ) µ(s
                                         0
                                         ) · l(s
                                              0
                                              ) // iterate lower vector
 9 u(s) := optµ∈T (s)
                             P
                               s0∈spt(µ) µ(s
                                         0
                                          ) · u(s
                                               0
                                                ) // iterate upper vector
10 return 1
             2
               (u(sI ) − l(sI ))
```
Alg. 1: Interval iteration for probabilistic reachability

#### 2.1 Algorithms

Interval iteration [3,5,12,13] computes reachability probabilities p(s) = P<sup>s</sup> opt(G), opt ∈ { max, min }. We show the basic algorithm as Alg. 1. It iteratively refines vectors l and u that map each state to a value in Q such that, at all times, we have l(s) ≤ p(s) ≤ u(s). In each iteration, the values in l and u are updated for all relevant states (line 7) via the classic Bellman equations of value iteration (lines 8-9). Their least fixpoint is p, towards which l converges from below. Some preprocessing is needed to ensure that the fixpoint is unique and also u converges towards p: for maximisation, we need to collapse MECs into single states (line 2). This can be be done via graph-based algorithms (see e.g. [7]) that only consider the graph structure of the MDP as in Definition 1 but do not perform calculations with the concrete probability values. For both maximisation and minimisation, we need to identify the sets S<sup>0</sup> and S<sup>1</sup> such that ∀s ∈ S<sup>0</sup> : p(s) = 0 and ∀s ∈ S<sup>1</sup> : p(s) = S<sup>1</sup> (line 3). This can equally done via graph-based algorithms [10, Algs. 1-4]. We then initialise l and u to trivial under-/overapproximations of p (lines 4-5). Iteration stops when the relative difference between l and u at s<sup>I</sup> is at most (which is often chosen as 10<sup>−</sup><sup>3</sup> or 10<sup>−</sup><sup>6</sup> ). The corresponding check in line 6 assumes that division by zero results in +∞, as is the default in IEEE 754. By convergence of l and u towards the fixpoint, II terminates, and we eventually return a value pˆ with the guarantee that p(s<sup>I</sup> ) ∈ [(1 − ) · p, ˆ (1 + ) · pˆ]. This makes II sound.

PCTL. The temporal logic PCTL [16] allows us to construct complex branchingtime properties. It takes standard CTL [2] and replaces the A(ψ) ("for all paths ψ holds") and E(ψ) ("there exists a path for which ψ holds") operators by the probabilistic operator P<sup>∼</sup>c(ψ) for "under all schedulers, the probability of the measurable set of paths for which ψ holds is ∼ c" where ∼ ∈ { <, ≤, >, ≥ } and c ∈ [0, 1]. To model-check a PCTL formula on MDP M, we follow the standard recursive CTL model checking algorithm [2, Sect. 6.4] except for the P