# Bridging Impulse Control of Piecewise Deterministic Markov Processes and Markov Decision Processes: Review and Open Challenges

Alice Cleynen, Benoîte de Saporta, Orlane Rossini, Régis Sabbadin, Amélie Vernay

► To cite this version:

Alice Cleynen, Benoîte de Saporta, Orlane Rossini, Régis Sabbadin, Amélie Vernay. Bridging Impulse Control of Piecewise Deterministic Markov Processes and Markov Decision Processes: Review and Open Challenges. 2026. hal-05458302

HAL Id: hal-05458302

https://hal.science/hal-05458302v1

Preprint submitted on 14 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Bridging Impulse Control of Piecewise Deterministic Markov Processes and Markov Decision Processes: Review and Open Challenges*

Alice Cleynen††

Benoîte de Saporta‡

Orlane Rossini‡

Régis Sabbadin§

Amélie Vernay‡

## Abstract

Control theory plays a pivotal role in understanding and optimizing the behavior of complex dynamical systems across various scientific and engineering disciplines. Two key frameworks that have emerged for modeling and solving control problems in stochastic systems are Piecewise Deterministic Markov Processes (PDMPs) and Markov Decision Processes (MDPs). Each framework has its unique strengths, and their intersection offers promising opportunities for tackling a broad class of problems, particularly in the context of impulse control and sequential decision-making for complex systems. The relationship between PDMPs and MDPs is a natural subject of exploration, as embedding impulse control problems for PDMPs into the MDP framework could open new avenues for their analysis and resolution. Specifically, this integration would allow leveraging the computational and theoretical tools developed for MDPs to address the challenges inherent in PDMPs. On the other hand, PDMPs can offer a versatile and simple paradigm to model continuous time problems that are often described as discrete-time MDPs parametrized by complex transition kernels. This transformation has the potential to bridge the gap between the two frameworks, enabling solutions to previously intractable problems and expanding the scope of both fields. This paper presents a comprehensive review of two research domains, illustrated through a running medical example. The example is revisited and progressively formalized within the framework of the various concepts and objects introduced.

## Contents

1 Introduction 3

2 Piecewise Deterministic Markov Processes 4

<table><tr><td>2.1</td><td>PDMP, definition and examples</td><td>4</td></tr><tr><td>2.1.1</td><td>Generic definition</td><td>4</td></tr><tr><td>2.1.2</td><td>Running medical example as a PDMP</td><td>6</td></tr><tr><td>2.1.3</td><td>Semi-Markov dynamics</td><td>7</td></tr><tr><td>2.1.4</td><td>Running medical example as a time-augmented PDMP</td><td>8</td></tr><tr><td>2.1.5</td><td>Simulation</td><td>8</td></tr><tr><td>2.1.6</td><td>Embedded chains</td><td>10</td></tr><tr><td>2.1.7</td><td>PDMPs as versatile models</td><td>10</td></tr><tr><td>2.1.8</td><td>Main research domains involving PDMPs</td><td>11</td></tr><tr><td>2.2</td><td>Impulse control for PDMPs</td><td>18</td></tr><tr><td>2.2.1</td><td>Continuous versus impulse control</td><td>18</td></tr><tr><td>2.2.2</td><td>Definition of an impulse control problem</td><td>18</td></tr><tr><td>2.2.3</td><td>Running medical example as a controlled PDMP</td><td>20</td></tr><tr><td>2.2.4</td><td>Resolution of impulse control problems for PDMP</td><td>21</td></tr><tr><td>2.2.5</td><td>Partially observed controlled PDMPs</td><td>22</td></tr></table>

<footnote>*This work was supported by the French Agence Nationale de la Recherche (ANR), grant ANR-21-CE40-005 (project HSMM-INCA) and the European Union's Horizon 2020 research and innovation program (Marie Sklodowska-Curie grant agreement No 890462)</footnote>

<footnote>†John Curtin School of Medical Research, Australian National University</footnote>

<footnote>‡IMAG, Univ. Montpellier, CNRS, Montpellier, France</footnote>

<footnote>§INRAE-MIAT, Univ Toulouse, Toulouse, France</footnote>
<table>
  <tbody>
    <tr>
      <td>3</td>
      <td>Markov Decision Processes</td>
      <td>22</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Markov Decision Process definition and control</td>
      <td>22</td>
    </tr>
    <tr>
      <td>3.1.1</td>
      <td>Generic definition</td>
      <td>23</td>
    </tr>
    <tr>
      <td>3.1.2</td>
      <td>MDP policies</td>
      <td>24</td>
    </tr>
    <tr>
      <td>3.1.3</td>
      <td>Policy Evaluation criteria</td>
      <td>24</td>
    </tr>
    <tr>
      <td>3.1.4</td>
      <td>Computing optimal policies</td>
      <td>25</td>
    </tr>
    <tr>
      <td>3.1.5</td>
      <td>Running medical example as a finite state and action MDP</td>
      <td>27</td>
    </tr>
    <tr>
      <td>3.1.6</td>
      <td>Generalizations</td>
      <td>28</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Partially Observed Markov Decision Processes</td>
      <td>29</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>POMDP definition</td>
      <td>29</td>
    </tr>
    <tr>
      <td>3.2.2</td>
      <td>Belief MDP</td>
      <td>31</td>
    </tr>
    <tr>
      <td>3.2.3</td>
      <td>Running medical example as a POMDP</td>
      <td>32</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Bayesian approach to handle MDPs with unknown model</td>
      <td>32</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Bayes-Adaptive MDP</td>
      <td>33</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Running medical example as a BAMDP</td>
      <td>34</td>
    </tr>
    <tr>
      <td>3.3.3</td>
      <td>Bayes-Adaptive POMDP</td>
      <td>35</td>
    </tr>
    <tr>
      <td>3.3.4</td>
      <td>Running medical example as a BAPOMDP</td>
      <td>36</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Algorithms for computing exact or approximate optimal policies in MDPs</td>
      <td>36</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Dynamic programming algorithms</td>
      <td>37</td>
    </tr>
    <tr>
      <td>4.1.1</td>
      <td>Backward iteration</td>
      <td>37</td>
    </tr>
    <tr>
      <td>4.1.2</td>
      <td>Value iteration</td>
      <td>38</td>
    </tr>
    <tr>
      <td>4.1.3</td>
      <td>Policy iteration</td>
      <td>38</td>
    </tr>
    <tr>
      <td>4.1.4</td>
      <td>Limitations and extensions</td>
      <td>38</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Monte-Carlo simulations-based algorithms</td>
      <td>39</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Reinforcement-learning algorithms</td>
      <td>41</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Model-free approaches: value-based algorithms</td>
      <td>42</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Model-free approaches: policy-based algorithms</td>
      <td>43</td>
    </tr>
    <tr>
      <td>4.3.3</td>
      <td>Model-based algorithms</td>
      <td>44</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Controlled PDMP as member of the MDP family</td>
      <td>46</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Controlled PDMP as MDP</td>
      <td>46</td>
    </tr>
    <tr>
      <td>5.1.1</td>
      <td>Impulse control strategies and MDP policies</td>
      <td>46</td>
    </tr>
    <tr>
      <td>5.1.2</td>
      <td>Formalism of the equivalent MDP</td>
      <td>47</td>
    </tr>
    <tr>
      <td>5.1.3</td>
      <td>Running medical example as a controlled PDMP formalized as an MDP</td>
      <td>49</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Partially observed controlled PDMP as POMDP</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>Formulation as a POMDP</td>
      <td>51</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>Running medical example as a partially observed controlled PDMP formalized as a POMDP</td>
      <td>52</td>
    </tr>
    <tr>
      <td>5.2.3</td>
      <td>Belief states for PDMP</td>
      <td>53</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Concluding remarks and open questions</td>
      <td>54</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Open questions in controlled PDMPs that might be tackled using MDPs</td>
      <td>54</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Interesting questions in MDPs arising from converted PDMPs</td>
      <td>55</td>
    </tr>
  </tbody>
</table>
# 1 Introduction

Control theory plays a pivotal role in understanding and optimizing the behavior of complex dynamical systems across various scientific and engineering disciplines. Two key frameworks that have emerged for modeling and solving control problems in stochastic systems are Piecewise Deterministic Markov Processes (PDMPs) and Markov Decision Processes (MDPs). Each framework has its unique strengths, and their intersection offers promising opportunities for tackling a broad class of problems, particularly in the context of impulse control and sequential decision-making for complex systems.

PDMPs, introduced in the 80s by Davis[44], are a versatile class of continuous-time and hybrid (discrete and continuous) state space stochastic processes that combine deterministic dynamics with discrete stochastic events, allowing for a realistic representation of systems with hybrid behavior. Impulse control problems for PDMPs are of particular interest in applications ranging from resource allocation and maintenance planning to biological systems and energy management. These problems involve determining the optimal timing and magnitude of interventions, often modeled as impulses, to minimize a cost or maximize a reward. However, the complexity of PDMPs, particularly their reliance on continuous state spaces and non-trivial event-driven dynamics, poses significant challenges for both theoretical analysis and numerical resolution of impulse control problems.

On the other hand, MDPs, introduced in the 50s by Bellman [14], provide a well-established framework for sequential decision-making under uncertainty, originally leveraging the discrete state and action spaces to facilitate analytical and computational approaches. They have been extended to very general state and action spaces. Extensions of the classical MDP framework, such as semi-Markov Decision Processes, Partially Observed MDPs or Bayesian MDPs offer increased flexibility in modeling various real-world problems. MDPs have been instrumental in advancing the study of stochastic control, particularly through their compatibility with dynamic programming and approximate inference techniques.

The relationship between PDMPs and MDPs is a natural subject of exploration, as embedding impulse control problems for PDMPs into the MDP framework could open new avenues for their analysis and resolution. Specifically, this integration would allow leveraging the computational and theoretical tools developed for MDPs to address the challenges inherent in PDMPs. On the other hand, PDMPs can offer a versatile and simple paradigm to model continuous time problems that are often described as discrete-time MDPs parametrized by complex transition kernels. This transformation has the potential to bridge the gap between the two frameworks, enabling solutions to previously intractable problems and expanding the scope of both fields.

This paper presents a comprehensive review of two research domains, illustrated through a running medical example. The example is revisited and progressively formalized within the framework of the various concepts and objects introduced. The paper is organized as follows.

* Section 2 provides an introduction to PDMPs and their impulse control problems, highlighting the mathematical structure and the challenges posed by their hybrid nature.

* Section 3 delves into the framework of MDPs and their extensions, introducing control policies and formulation of MDPs under partial observation or unknown models.

* Section 4 presents the state-of-the-art in solving decision-making problems under uncertainty by introducing tools and algorithms relevant to control problems formulated as MDPs.

* Section 5 demonstrates how impulse control problems for PDMPs can be embedded into the MDP framework. We also identify a series of open questions that could be addressed through this integration, such as scalability, approximate solutions, and new modeling paradigms.

By bringing these two frameworks together, we aim to contribute to the growing field of research at the interface of stochastic control, decision theory, statistics, computational sciences and reinforcement learning, offering insights into how embedding PDMPs into the MDP framework can advance the field and tackle outstanding challenges in control and optimization. Throughout this paper, we use a medical case study as a guiding example to illustrate the theoretical concepts and methodologies being introduced. This red thread example will evolve progressively in complexity as new notions are presented. Starting with a simple scenario involving the free dynamics of a cancer marker, we will incrementally incorporate more realistic elements such as treatment or surgery impact on the marker value, and uncertain observations. For the simple scenarios, a Jupyter notebook is provided with this survey to illustrate simulations and resolutions. By building on the same context throughout the paper, we aim to provide readers with a
consistent and practical framework to understand and apply the concepts discussed, ultimately showcasing
the power of the proposed approaches in addressing real-world challenges.

## 2 Piecewise Deterministic Markov Processes

Piecewise Deterministic Markov Processes (PDMPs) are a general class of non-diffusion processes introduced by M. Davis in the 80's (see [44]) covering a wide range of applications, from workshop optimization, queuing theory, internet networks, reliability, insurance and finance or biology for instance. These continuous-time processes extend the class of point processes [45, 85, 35] and are characterized by deterministic motions punctuated by random jumps. They are especially suitable to model complex systems with interacting random phenomena, or large population asymptotics.

This presentation of PDMPs is divided into two main topics. In Section 2.1 we give the definition and examples of application of PDMPs, and provide their key features such as simulation properties. In Section 2.2, we provide a survey on the impulse control problem for this class of processes.

### 2.1 PDMP, definition and examples

We start with the definition and generic properties of PDMPs. In Section 2.1.1, we provide a concise formal definition of PDMPs. We introduce our running example, which will serve as a reference to illustrate the models discussed throughout the paper, in Section 2.1.2. Extensions to semi-Markov dynamics are explored in Section 2.1.3, with further application to the running example detailed in Section 2.1.4. We then delve into key questions and tools related to PDMPs, beginning with simulation techniques in Section 2.1.5. The identification and analysis of specific sub-chains of interest are addressed in Section 2.1.6. Representative examples of PDMPs from their main application domains are presented in Section 2.1.7, followed by an overview of prominent research questions and challenges in Section 2.1.8.

#### 2.1.1 Generic definition

We start with a formal definition of a PDMP, its state space and its local characteristics.

**Definition 1 (PDMP).** A Piecewise Deterministic Markov Process (PDMP) $X = (X_t)_{t \ge 0}$ is defined by a tuple $(E, \Phi, \lambda, Q)$, where

* The state space $E$ is a finite union of Borel subsets of $\mathbb{R}^d$, for some $d \ge 1$ and $\mathcal{B}(E)$ is its Borel $\sigma$-field. It is the set of all possible states of the process. Denote $\bar{E}$ the closure of $E$ and $\mathcal{B}(\bar{E})$ the Borel $\sigma$-field on $\bar{E}$.

* The flow $\Phi$ is a continuous function from $E \times \mathbb{R}_+$ onto $\bar{E}$ satisfying a semi-group property, i.e. $\Phi(\cdot, t + s) = \Phi(\Phi(\cdot, t), s)$ for all $s, t \in \mathbb{R}_+$. The flow prescribes the deterministic motion between jumps. It may be defined by an explicit function or simply described as the solution to a system of ordinary differential equations (ODEs).

* The jump intensity, or hazard rate, or risk function $\lambda$ is a measurable function from $\bar{E}$ onto $\mathbb{R}_+$ that determines the occurrence of random jumps.

* The jump kernel $Q$ is a Markov kernel on $(\mathcal{B}(E), \bar{E})$ that selects the new location of the process after each jump.

The PDMP dynamics can be described informally as follows: starting from some initial point $x \in E$, the motion of the process follows the deterministic flow $t \mapsto \Phi(x, t)$ until a first jump time $T_1$, which may be deterministically triggered by the process reaching the boundary $\partial E$ of $E$, or may occur randomly with hazard rate $t \mapsto \lambda(\Phi(x, t))$. At $T_1$ the process jumps to a new point $x' = X_{T_1}$ selected with probability $Q(dx'|\Phi(x, T_1))$ (conditional on $T_1$), and the motion restarts from this new point as before. The flow $\Phi$, the jump rate $\lambda$, and the Markov kernel $Q$ are called the local characteristics of the PDMP.

A generic representation of a PDMP is given in Figure 1. It is particularly convenient, when dealing with PDMPs, to separate the state space $E$ and states $x$ into discrete and Euclidean components. In
<figure><img src="image_2.png" /><figcaption>Figure 1: **Trajectory a generic PDMP.** Starting from an initial value $x_0$ at time $t = 0$, the process follows a deterministic trajectory until a jump occurs, either at a random time (as at $T_1$) or because the process reaches the state boundary (as at $T_2$). At jump times, the process jumps to a new location drawn from kernel $Q$.</figcaption></figure>

this case the state space is called *hybrid*. Let $M$ be a finite set, and for all $m \in M$, let $E_m$ be some Borel subset of $\mathbb{R}^{d_m}$, which dimension $d_m$ is allowed to change with $m$. The state space is thus

$$E = \bigcup_{m \in M} \{m\} \times E_m.$$

Any state $x \in E$ can be written as $x = (m, x)$ with $m \in M$ and $x \in E_m$. The discrete component $m$ is called the *mode* or *regime* and $x$ is called the *Euclidean variable*. In between jumps, the mode $m$ remains constant and only the Euclidean variable $x$ evolves through the flow. At jump times, both the mode and the Euclidean variable may be allowed to change through the Markov kernel $Q$.

In the hybrid setting, we set

$$\Phi(x, t) = (m, \Phi_m(x, t)),$$

for all $x = (m, x) \in E$ and $t \in \mathbb{R}_+$. The jump intensity can be written as $\lambda(x) = \lambda_m(x)$ for all $x = (m, x) \in E$ and $t \in \mathbb{R}_+$ and is such that for any $x$ in $E$, there exists $\epsilon > 0$ such that

$$\int_{0}^{\epsilon} \lambda(\Phi(x, t))dt < +\infty,$$

forbidding instantaneous jumps. Let $\partial E$ be the boundary of $E$: $\partial E = \cup_{m \in M} \{m\} \times \partial E_m$, and $t^*(x)$ be the deterministic time the flow takes to reach the boundary $\partial E$ of $E$ when it starts from $x$

$$t^*(x) = t_m^*(x) = \inf\{t > 0 : \Phi(x, t) \in \partial E\}.$$

Jumps may occur via two means: *random jumps* occur from the realization of the random clock with intensity $\lambda$ (e.g. time $T_1$ on Figure 1), while *deterministic jumps* occur when the process reaches the boundary of the state space (e.g. time $T_2$ on Figure 1). Thus, starting from $x$ at time 0, the first jump time $T_1$ has the following distribution

$$P_x(T_1 > t) = P(T_1 > t | X_0 = x) = \exp\left(-\int_0^t \lambda(\Phi(x, s))ds\right) 1_{t < t^*(x)}, \quad t \ge 0.$$

For convenience, we also define

$$\Lambda(x, t) = \Lambda_m(x, t) = \int_0^t \lambda(\Phi(x, s))ds,$$
**Algorithm 1** Iterative construction of a trajectory of a PDMP starting from state $x_0 = (m_0, x_0)$ at time 0 and up to the N-th jump time

1: $x \leftarrow x_0$

2: $T \leftarrow 0$

3: **for** $n \leftarrow 0$ to $N-1$ **do**

4: $S \sim \lambda(x)$

5: $S \leftarrow \min\{S, t^*(x)\}$

6: $X_t \leftarrow \Phi(x, t-T)$ for $T \le t < T+S$

7: $x' \sim Q(\cdot|\Phi(x, S))$

8: $x \leftarrow x'$

9: $T \leftarrow T+S$

10: **end for**

the cumulative risk up to $t$ starting from $x$. Finally, the Markov kernel $Q$ on $(\mathcal{B}(E), \bar{E})$ selects the new location after each jump. We also set

$$Q_m(\cdot|x) = Q(\cdot|x),$$

for $x = (m, x)$ in accordance with other notations. It satisfies $Q(\{x\}|x) = 0$ so that a jump has to change the state of the process (either the mode, or the Euclidean variable or both).

After the first jump time $T_1$, the process restarts from a new state selected by $Q$ and follows the flow until the next jump time $T_2$. The sojourn time $T_2 - T_1$ until the next jump is drawn from the jump intensity and time to reach the boundary as before. At the jump time $T_2$, the post-jump location is selected by the Markov kernel $Q$ and so on iteratively, as described in Algorithm 1. Note that Algorithm 1 can only simulate the process up to time $T_\infty = \lim_{n \to \infty} T_n$ that might be finite. This phenomenon is called explosion. To avoid explosion, one usually also requires that for all $x \in E$, one has

$$\mathbb{E}_{x}[\lim_{n \to \infty} T_n] = +\infty.$$

A process $X$ defined by a flow $\Phi$, an intensity function $\lambda$ and a kernel $Q$ verifying the above properties is well defined and is strong Markov, as is shown, for instance, in [45] or [85].

### 2.1.2 Running medical example as a PDMP

Let us now give a simple example of PDMP. This example is a guiding thread. It will be modified, adapted and enriched to illustrate the main families of models exposed throughout this paper.

We consider a patient undergoing medical follow-up after developing cancer. A biomarker level is measured during medical visits as a proxy for the progression of the disease. The evolution of the biomarker is modeled through a PDMP $X_t = (m_t, \zeta_t)$, where $m$ will denote the condition of the patient and $\zeta$ the marker level. At the start of monitoring, the patient is administered a treatment, the effect of which is to reduce the biomarker level exponentially fast. If the level reaches a certain fixed threshold $\zeta_0$, the patient is considered to be in remission and the marker remains constant. In the (random) event of a relapse, the biomarker level rises again exponentially fast. A sample trajectory is illustrated in Figure 2.

Thus, there are three possible modes for the patient corresponding to their condition. They can either be sick under treatment ($m = -1$), sick without treatment ($m = 1$) or in remission ($m = 0$). For simplicity, we assume that the process characteristics in remission mode are the same with or without treatment and do not differentiate the two cases. Set $M = \{-1, 0, 1\}$ and $E_{-1} = (\zeta_0, +\infty)$, $E_0 = \{\zeta_0\}$ and $E_1 = (\zeta_0, +\infty)$. The state space of this PDMP is

$$E = \bigcup_{m \in M} \{m\} \times E_m.$$

It is a hybrid state space with $M = M$, $m = m$, $x = \zeta$ and $E_m = E_m$ for all $m$.

* Under treatment $m = -1$, the biomarker level $\zeta$ decreases exponentially with a slope $v_{-1} < 0$. During relapse $m = 1$, it increases exponentially with a slope $v_1 > 0$. For all $\zeta \in \mathbb{R}$ and for all
<figure><img src="image_3.png" /><figcaption>Figure 2: **Trajectory of the running medical example as a PDMP.** The PDMP starts in mode $m = -1$ (patient under treatment) from an initial point $x_0$ and follows a deterministic trajectory along its exponential flow $\Phi_{-1}$ until the first jump occurs at time $T_1$ when reaching the boundary $x = \zeta_0$ (remission). The mode switches to $m = 0$ and the flow is constant equal to $\zeta_0$ until a new jump occurs at time $T_2$ with an exponential clock. The mode switches then to $m = 1$ (relapse of the disease) and the trajectory rises exponentially along the flow $\Phi_1$.</figcaption></figure>

$t \in \mathbb{R}_+$ we have

$$ \left\{ \begin{array}{l} \Phi_{-1}(\zeta, t) = \zeta e^{v_{-1}t}, \\ \Phi_0(\zeta, t) = \zeta = \zeta_0, \\ \Phi_1(\zeta, t) = \zeta e^{v_1 t}. \end{array} \right. $$

* In mode $m = -1$, a jump occurs when the subject reaches the fixed remission threshold $\zeta_0$. This is a deterministic jump at the boundary and $\lambda_{-1}(\zeta) = 0$. The jump time $t^*_{-1}(\zeta)$ is the solution of $\Phi_{-1}(\zeta, t) = \zeta_0$. That is, $t^*_{-1}(\zeta) = \frac{1}{v_{-1}}\log(\frac{\zeta_0}{\zeta})$. In mode $m = 0$, jumps may only occur randomly, $t^*_0(\zeta) = +\infty$ and the jump intensity $\lambda_0$ is a positive constant. Once the process reaches mode $m = 1$, no more jump can occur and $\lambda_1(\zeta) = 0$, $t^*_1(\zeta) = +\infty$.

* The mode-specific Markov kernels are

$$ \left\{ \begin{array}{l} Q(m', \zeta' | (-1, \zeta)) = 1_{\zeta=\zeta_0=\zeta'} \times 1_{m'=0}, \\ Q(m', \zeta' | (0, \zeta)) = 1_{\zeta=\zeta_0=\zeta'} \times 1_{m'=1}. \end{array} \right. $$

### 2.1.3 Semi-Markov dynamics

In the example above, when the patient is in mode 0, the flow is constant. Hence, the jump intensity also has to be constant, i.e. it corresponds to a memory-less exponential distribution. In this case, as in many other practical examples, see e.g. [11], the exponential distribution may not be realistic to model the underlying phenomenon. However allowing for more general distributions breaks the Markov property and falls into the class of semi-Markov processes.

It is especially easy to encompass semi-Markov dynamics into the PDMP framework by enlarging the state space. This new process is sometimes called the *time-augmented process*, see [45, Section 31.5] for a slightly different definition in the same spirit.

**Definition 2 (Time-augmented PDMP).** A time-augmented Piecewise Deterministic Markov Process is a PDMP defined by the tuple $(E, \Phi, \lambda, Q)$, where
* The state space E has the specific form

$$E = \bigcup_{m \in M} \{m\} \times E_m \times \mathbb{R}_+,$$

for some finite mode set M, and where for all $m \in M$, $E_m$ is some Borel subset of $\mathbb{R}^{d_m}$.

* The flow $\Phi$ has the specific form

$$\Phi(x, t) = (m, \Phi_m(x), u + t),$$

for $x = (m, x, u) \in E$, where $\Phi_m$ are flows from $E_m \times \mathbb{R}_+$ onto $E_m$ satisfying the semi-group property.

* The jump intensity just satisfies the conditions from definition 1 and is allowed to depend on all coordinates of $x \in E$.

* The jump kernel $Q$ has the specific form

$$Q(\{m'\} \times A \times B|x) = \mathbb{1}_{0 \in B} Q_m(A|x, u),$$

for $x = (m, x, u) \in E$, A any Borel subset of $\cup\{m\} \times E_m$, B any Borel subset of $\mathbb{R}_+$, and where $Q_m$ is a Markov kernel on $(B(\cup\{m\} \times E_m), \bar{E}_m)$.

The additional component $u$ can be interpreted as the time elapsed since the last jump. It moves along the flow at speed 1 and is reset to 0 at every jump.

One could define more sophisticated time-augmented processes where either $u$ is never reset to 0 (it is thus the time since the beginning, as in [45, section 31.5]), resets less frequently than at each jump, or is multidimensional, see Section 2.1.7 for such examples.

### 2.1.4 Running medical example as a time-augmented PDMP

We consider the same medical example as in Section 2.1.2, but in which we want to allow the risk of relapse to depend on how long the patient has been in the remission phase. To this effect, we introduce the extra variable $u$ in mode $m = 0$ representing the time elapsed since the entry into remission (i.e. since the last jump to mode $m = 0$). In mode 0, the state space thus becomes

$$E_0 = \{\zeta_0\} \times \mathbb{R}_+.$$

In modes $m = -1$ and $m = 1$, the state spaces, flows and jump intensities remain unchanged and the mode-specific Markov kernel becomes $Q(m', \zeta', u' | (-1, \zeta)) = \mathbb{1}_{\zeta=\zeta_0=\zeta'} \times \mathbb{1}_{u'=0} \times \mathbb{1}_{m'=0}$.

In mode $m = 0$, the flow is defined by $\Phi_0(\zeta, u, t) = (\zeta_0, u + t)$, the intensity may now depend on $u$, for instance through a Weibull distribution i.e. $\lambda_0(\zeta, u) = \beta u^\alpha$, and the Markov kernel becomes $Q(m', \zeta' | (0, \zeta, 0)) = \mathbb{1}_{\zeta=\zeta_0=\zeta'} \times \mathbb{1}_{m'=1}$.

### 2.1.5 Simulation

A very interesting property of PDMPs is that they are particularly easy to simulate due to their intrinsic iterative construction as seen in Algorithm 1.

Computing the trajectory along the flow (Alg. 1, line 6) is straightforward if the flow is explicit. If it is only known as the solution to a system of ODEs, one just needs to use any suitable ODE solver.

Computing the post-jump location from the kernel $Q$ (Alg. 1, line 7) is also usually straightforward. In most applications, the kernel only changes the mode and leaves the Euclidean variable unchanged, or changes the Euclidean variable through some known mapping depending only on its value just before the jump, see examples in Section 2.1.7.

The only challenging step is then to simulate the jump dates from the non constant intensity (Alg. 1, line 4). We present here an exact simulation method called the Stochastic Simulation Algorithm (SSA), or Gillespie algorithm [69] that is especially suitable for PDMPs. The interested reader may consult [34, 46] for an alternative algorithm to simulate only sojourn-times and post-jump locations instead of full trajectories, and [108, 127] and references therein for a broader discussion of alternative exact and approximate simulation methods.

The SSA is based on the Markov property and the following property of the minimum of independent generalized exponential distributions.
**Algorithm 2** Stochastic Simulation Algorithm for a PDMP starting from state $x_0 = (\mathbf{m}_0, \mathbf{x}_0)$ at time 0 and up to time horizon $H$

1: $x \leftarrow x_0$

2: $t \leftarrow 0$

3: **while** $t < H$ **do**

4: $S \sim \text{Exp}(\bar{\lambda})$

5: $X_{t+s} \leftarrow \Phi(x, s)$ for $0 \le s < \min\{S, t^*(x)\}$

6: **if** $S > t^*(x)$ **then**

7: $x' \sim Q(\cdot|\Phi(x, t^*(x)))$

8: **else**

9: $U \sim U[0, 1]$

10: **if** $U \le \frac{\lambda(S)}{\bar{\lambda}}$ **then**

11: $x' \sim Q(\cdot|\Phi(x, S))$

12: **else**

13: $x' \leftarrow \Phi(x, S)$

14: **end if**

15: **end if**

16: $x \leftarrow x'$

17: $t \leftarrow t + \min\{S, t^*(x)\}$

18: **end while**

**Lemma 1.** Let $S_1$ and $S_2$ be two independent random variables with respective (non-constant) intensity functions $\lambda_1$ and $\lambda_2$. Set $S = \min\{S_1, S_2\}$ and $I = \arg\min\{S_1, S_2\}$. Then the intensity of $S$ is $\lambda_1 + \lambda_2$ and conditionally to $(S = t)$, the distribution of $I$ is $(\frac{\lambda_1(t)}{\lambda_1(t)+\lambda_1(t)}, \frac{\lambda_2(t)}{\lambda_1(t)+\lambda_1(t)})$.

*Proof.* First, notice that if $S_i$ has intensity $\lambda_i$, then its density function is $t \mapsto \lambda_i(t)e^{-\int_0^t \lambda_i(s)ds}$. Using the independence between $S_1$ and $S_2$, one obtains the joint distribution of $(S, I)$ as follows

$$ 
\begin{aligned}
\mathbb{P}(S > t, I = 1) &= \mathbb{P}(S_1 > t, S_2 > S_1) \\
&= \int_t^\infty \lambda_1(x) e^{-\int_0^x \lambda_1(s)ds} \int_x^\infty \lambda_2(y) e^{-\int_0^y \lambda_2(s)ds} dy dx \\
&= \int_t^\infty \lambda_1(x) e^{-\int_0^x \lambda_1(s)ds} e^{-\int_0^x \lambda_2(s)ds} dx \\
&= \int_t^\infty \frac{\lambda_1(x)}{\lambda_1(x) + \lambda_2(x)} (\lambda_1(x) + \lambda_2(x)) e^{-\int_0^x (\lambda_1(s)+\lambda_2(s))ds} dx.
\end{aligned}
$$

Hence the density of $S$ is $(\lambda_1(x) + \lambda_2(x))e^{-\int_0^x (\lambda_1(s)+\lambda_2(s))ds}$, and

$$ \mathbb{P}(I = 1 | S = t) = \frac{\lambda_1(x)}{\lambda_1(x) + \lambda_2(x)}, $$

and the result follows.

<figure><img src="image_4.png" /></figure>

The principle underlying the SSA algorithm is as follows. The aim is to sample from a random variable $S_1$ with (non-constant) intensity $\lambda_1$. Suppose that $\lambda_1$ is bounded from above by some constant $\bar{\lambda}$ and let $S$ be random variable with (regular) exponential distribution with parameter $\bar{\lambda}$, so that it is very easy to sample from $S$. In view of Lemma 1, $S$ can be interpreted as the minimum of $S_1$ and some virtual extra random variable $S_2$ with intensity $\lambda_2 = \bar{\lambda} - \lambda_1$. Let $t$ be a sample value from $S$, then $t$ can be considered as a value sampled from $S_1$ with probability $\frac{\lambda_1}{\bar{\lambda}}$. This simple fact leads to a rejection algorithm where one repeatedly samples from $S$ until the result is accepted as a sample from $S_1$. In a Markov context, it is even more efficient as by construction $S_1 \ge S$. Hence all samples are used, some are accepted as jump times (boundary jump: Alg. 2 line 6 or random jump: Alg. 2 line 10), the others just move the process forward without creating a jump (Alg. 2 line 12). The procedure then restarts from the new location, as detailed in Algorithm 2.

The choice of the upper bound $\bar{\lambda}$ can have a strong impact on the speed of the algorithm. If it is too high, it will keep proposing very short times and move forward extremely slowly. If at all possible, one should choose a tight upper bound for the intensity.
### 2.1.6 Embedded chains

As for any continuous-time Markov process, one can define many different discrete-time Markov chains embedded in a PDMP. Some of them are particularly interesting either for theory or applications and are specified below.

**Canonical chain.** The canonical chain is the most natural Markov chain embedded in a PDMP. It corresponds to the post-jump-locations and inter-jump-times chains. Let us define $Z_0 = X_0$ and $S_0 = 0$, and then iteratively, for all $n \in \mathbb{N}^*$, $Z_n = X_{T_n}$ the value of the PDMP at the $n^{th}$ jump time and $S_n = T_n - T_{n-1}$ the inter-jump time between the $n-1^{th}$ and $n^{th}$ jumps. Then $(Z_n, S_n)_{n \in \mathbb{N}}$ is a Markov chain called the *canonical embedded chain* associated to the PDMP process $X$. Its transition kernel $P^{Z,S}$ is defined, for any Borel subsets $A$ of $E$ and $B$ of $\mathbb{R}_+$, by

$$ P^{Z,S}(A \times B|x, t) = P^{Z,S}(A \times B|x) \\ = \int_{0}^{t^*(x)} \lambda(\Phi(x, s))e^{-\Lambda(x, s)}Q(A|\Phi(x, s))1_{B}(s)ds + e^{-\Lambda(x, t^*(x))}Q(A|\Phi(x, t^*(x)))1_{B}(t^*(x)), $$

where the first part of the equation corresponds to random jumps, while the second part corresponds to deterministic jumps.

As the motion of the PDMP is deterministic between jumps, it is straightforward to reconstruct the full trajectory of the continuous-time PDMP from the knowledge of $(Z_n, S_n)_{n \in \mathbb{N}}$ as

$$ X_t = \Phi \left( Z_n, t - \sum_{i=0}^{n} S_i \right) \quad \text{if } t \in \left[ T_n = \sum_{i=0}^{n} S_i, T_{n+1} = \sum_{i=0}^{n+1} S_i \right), $$

with the convention $T_0 = S_0 = 0$. Thus, the canonical chain contains all the information of the process and is especially useful to obtain recursive formulations based on the Markov property. It is also at the heart of several numerical approximations for PDMPs, see e.g. [46].

For some applications, especially in inference and MCMC frameworks (see Section 2.1.8), this canonical chain can be extended with the pre-jump locations $\Phi(Z_{n-1}, S_n)$. In some numerical procedures, it can also be extended with the jump times $T_n$. All these chains are especially interesting as they capture all the randomness of the process and have explicit transition kernels with a relatively simple form depending on the local characteristics of the PDMP.

**Skeleton chains.** Other common embedded chains are the *skeleton chains* defined by sampling the PDMP on a deterministic (but not necessarily homogeneous) grid $(t_1, \dots, t_n)$. We simply denote $X_n = X_{t_n}$, which once again is a (time inhomogeneous) Markov chain with kernels $P_{t_{n+1}-t_n}$.

While kernels $P_{t_{n+1}-t_n}$ might still have explicit forms, unless constraints are added on the PDMP there is no a priori limit in the number of jumps occurring between two time-points $t_n$ and $t_{n+1}$, which prevents a generic formulation of the kernels. Moreover, unless jump-points are observed, it is not possible to reconstruct the original PDMP process from the observation of $X_n$. However, in many application settings, $X_n$ might be the only available information from the process, see for instance the version of the running medical example developed in Section 5.2.2.

**Augmented canonical chain.** Interestingly, in [40], to improve comparisons of the canonical chain to the original process when jumps are scarce, the authors combine the canonical chain with a skeleton chain where the grid is randomly chosen. In this setting, the observations of the skeleton chain are constrained to occur in between jump times on a time-frame defined by an independent Poisson process of rate 1. Among other properties, this allows to translate necessary and sufficient conditions on the PDMP regarding reducibility, recurrence, existence and convergence to invariant measures to the resulting discrete-time Markov chain for which assumptions are much easier to verify.

### 2.1.7 PDMPs as versatile models

PDMPs are constructed using a small set of intuitive and easily interpretable local characteristics, making them a highly versatile and powerful modeling tool. They extend the framework of discrete-time and discrete-state space Markov or semi-Markov chains, enabling the representation of a broad spectrum of
<table><caption>Table 1: Some examples of PDMP and their main features</caption><thead><tr><th rowspan="2">Figure</th><th>mode</th><th>Eucl. var.</th><th rowspan="2">boundaries</th><th rowspan="2">flow</th><th rowspan="2">time augmented</th></tr><tr><th>cardinal of M</th><th>dimension of E<sub>m</sub></th></tr></thead><tbody><tr><td>Fig. 3</td><td>2</td><td>-</td><td>no</td><td>constant</td><td>no</td></tr><tr><td>Fig. 4</td><td>(N + 1)<sup>2</sup></td><td>-</td><td>no</td><td>constant</td><td>no</td></tr><tr><td>Fig. 5</td><td>-</td><td>1</td><td>no</td><td>linear</td><td>no</td></tr><tr><td>Fig. 6</td><td>-</td><td>1</td><td>no</td><td>exponential</td><td>no</td></tr><tr><td>Fig. 7</td><td>2</td><td>1</td><td>no</td><td>linear</td><td>no</td></tr><tr><td>Fig. 8</td><td>5</td><td>1</td><td>no</td><td>linear</td><td>yes</td></tr><tr><td>Fig. 9</td><td>2</td><td>2</td><td>no</td><td>analytical expr.</td><td>no</td></tr><tr><td>Fig. 10</td><td>6</td><td>3</td><td>no</td><td>analytical expr.</td><td>yes</td></tr><tr><td>Fig. 11</td><td>N + 1</td><td>1</td><td>no</td><td>sol. of an EDO</td><td>no</td></tr><tr><td>Fig. 12</td><td>3</td><td>1</td><td>yes</td><td>exponential</td><td>no</td></tr></tbody></table>

complex phenomena across various domains.

**Sub-classes of PDMPs** PDMPs encompass and generalize many classes of Markov processes, including

* discrete-time and space Markov chains (corresponding to a time-augmented PDMP with no Euclidean variables, null intensity, jumps triggered at boundaries at each time unit);

* continuous-time and discrete state space Markov processes (corresponding to a PDMP with no Euclidean variables and constant intensity for each mode);

* semi-Markov discrete state space processes (corresponding to a time-augmented PDMP with no Euclidean variables and time-dependent intensity).

**Examples of phenomena modeled by PDMPs** Examples from the PDMP literature cover a wide range of application domains, from workshop optimization and queuing theory (e.g. [45]), internet networks (e.g. [12]), reliability (e.g. [46]), insurance and finance (e.g. [23]), neurosciences (e.g. [107, 109]), population dynamics (e.g. [60, 36]), genomics (e.g. [79]), or medical treatment (e.g. [30, 103]) to name just a few.

Here we provide a list of selected examples that cover several characteristics of PDMPs: presence or absence of mode and Euclidean variables, boundaries, type of flows, etc. Figures 3 to 12 illustrate some sampled trajectories from the corresponding PDMP models and Table 1 summarizes their main features. While far from being exhaustive, the list of examples has been chosen to emphasize the versatility of the PDMP framework:

* PDMPs with no Euclidean variable: Fig. 3 and 4,

* PDMPs with no mode variable: Fig. 5 and 6,

* PDMPs with mode and Euclidean variables but no boundary: Fig. 7 to 11,

* PDMPs with mode and Euclidean variables and boundaries: Fig. 12.

Some authors have also extended the definition of PDMPs to infinite-dimensional states, such as measure spaces, to study population dynamics with a varying number of individuals, see e.g. [75, 32]. Those types of PDMPs will not be discussed further.

### 2.1.8 Main research domains involving PDMPs

PDMPs are being studied by various communities in different frameworks. Among the main aspects present in the literature, we very briefly present here some references to existing surveys, and will especially focus on control of PDMPs in the sequel of this work, see Sections 2.2 and 5.
<figure><img src="image_5.png" /><figcaption>Figure 3: **Sample trajectory of the telegraph model.** The PDMP has 2 modes (0 and 1), no Euclidean variable, and no boundary. It switches from one mode to the other with a constant intensity depending on the mode.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: **Sample trajectory of the stochastic SIR model.** The PDMP has $(N+1)^2$ modes of the form $\mathbf{m} = (s, i)$ where $N$ is the total population size, $s$ the number of susceptible individuals, $i$ the number of infected individuals, and $r = N-s-i$ the number of removed individuals. The jump intensities are constant and depend on the mode. The only possible transitions correspond to the infection of a susceptible individual $\mathbf{m}' = (s', i') = (s-1, i+1)$ ($r$ is unchanged) or the removal of an infected individual $\mathbf{m}' = (s', i') = (s, i-1)$ ($r' = r+1$). There is no Euclidean variable and no boundary. The right-hand-side figure is a detailed section of the main graph on the left emphasizing piecewise constant trajectories.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: **Sample trajectory of the transmission control protocol (TCP) model.** The PDMP has no mode, a 1-dimensional Euclidean variable with linear flow corresponding to the available window size of a communications network, and no boundary. The jump intensity increases with the window size, and jumps correspond to random congestion resulting in a punctual decrease of the window size. The model is detailed in [25].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: **Sample trajectory of the individual cell growth and fragmentation model.** The PDMP has no mode, a 1-dimensional Euclidean variable with exponential flow corresponding to the size of a cell, and no boundary. The jump intensity increases with the cell size, and jumps correspond to divisions of the cell. The model is detailed in [60]</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: **Sample trajectory of the Zig-zag model.** The PDMP has 2 modes, a 1-dimensional Euclidean variable with linear flow and no boundary. The Euclidean variable corresponds to the distance between a bacterium and its nutrient, and the mode is its direction. The bacterium direction changes randomly with an intensity that depends whether the bacterium moves towards the nutrient (lower jump rate) or away from it (higher jump rate). The model is detailed in [67].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: **Sample trajectory of the optronic equipment model.** The PDMP has 5 modes, a 1-dimensional Euclidean variable with linear flow corresponding to a time-augmentation and no boundary. The modes correspond to different possible states of the equipment, with a nominal state, a degraded state and 3 possible different failure states. The Euclidean variable is the age of the equipment, and jump rates are age dependent. The model is detailed in [13]</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: **Sample trajectory of the gene expression model.** The PDMP has 2 modes, a 2-dimensional Euclidean variable with analytical flows, and no boundary. The mode is the state of the promoter (0:off or 1:on), the first Euclidean variable is the quantity of transcripted messenger RNA (it increases when the promoter is on, and decreases otherwise), and the second Euclidean variable is the quantity of translated protein (that depends on the quantity of mRNA). The promoter switches mode with constant mode-dependent intensities. The model is detailed in [79].</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: **Sample trajectory of the corrosion model.** The PDMP has 6 modes, a 3-dimensional Euclidean variable with analytical flows, no boundary but is time-augmented in some modes. Modes $m = (a, m)$ correspond to the presence ($a = 1$)/absence ($a = 0$) of an anti-corrosion protection and three possible environments $m \in \{1, 2, 3\}$. The succession of environments follows a deterministic cycle with random durations. The Euclidean variables are the thickness loss due to corrosion, the aggressiveness of the environment and the time spent with active anti-corrosion protection if $a = 1$. The model is detailed in [48]</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: **Sample trajectory of the Morris Lecar model.** The PDMP has $N + 1$ modes, a 1-dimensional Euclidean variable solution of an ODE, and no boundary. The mode corresponds to the number of open ionic channels along a neuron (for a total of $N$ channels), and the Euclidean variable is the electric potential of the neuron. The model is detailed in [94].</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: **Sample trajectory of the crack propagation model.** The PDMP has 3 modes and a 1-dimensional Euclidean variable exponential flow with a boundary. The modes correspond to material characteristics influencing the propagation of cracks, and the Euclidean variable is the crack size. The process stops when the crack size reaches a failure threshold. The model is detailed in [29].</figcaption></figure>

**Long-time behavior.** One of the first standard questions when studying stochastic processes is to understand their long time behavior, and establish conditions to ensure their convergence to some invariant or ergodic distribution. In [40], the authors establish necessary and sufficient conditions on existence of and convergence to ergodic distributions for PDMP based on properties of an augmented canonical discrete-time chain embedded in the PDMP, allowing to use the powerful theory of convergence to equilibrium for discrete-time Markov chains.

Subsequent works have focused on obtaining an explicit form for the ergodic distribution, and investigating the speed of convergence to equilibrium, see e.g. [64, 15, 33].

**Large population asymptotics.** Another branch of PDMPs is their use in modeling population dynamics, especially when several populations with different time scales are interacting. PDMPs can be used to model individual-based measure-valued populations or can be obtained as limiting processes in the large population scale, see e.g. [107, 75, 66, 33, 32, 105].

**Monte Carlo methods.** In the detailed survey [63], Fearnhead and co-authors exhibit PDMPs as the unifying framework of continuous-time Monte Carlo Markov chains (MCMC) and sequential Monte Carlo (SMC) methods, which in fact both involve simulating PDMPs. They appear as the natural building blocks of continuous-time Monte Carlo methods thanks to their low computational simulation properties (see Section 2.1.5) as only a finite number of events and transitions occur in any fixed time-interval. An estimate of the target distribution $\pi$ against any given function $g$ can easily be computed as the average of $g$ over the path of the PDMP, provided the PDMP has invariant distribution $\pi$. In this context, PDMPs have no modeling purpose, and the reader is left to [63] and references therein for more properties and results. More recent considerations may also be found in [17] or [102].

**Inference.** Statistical estimation for PDMPs has been studied in the literature, but remains a challenging task in practice. A general overview of recent methods for statistical parameter estimation in PDMPs is given in [5]. It is based on a wide range of application examples, in both parametric and non-parametric frameworks. While it is usually easy to estimate the flow of a PDMP, estimating the jump intensity or the kernel function can often be more difficult. We present here a few references dealing with the specific estimation of these characteristics in different contexts.

Examples of non parametric methods to estimate the jump mechanism under different settings and assumptions can be found in [6], [8] and [91]. In [6], the authors propose an estimation procedure for the conditional distribution of the inter-arrival times of a PDMP. In [91], the author provide a method to estimate the jump rate of a process which state space is (0, $\infty$) and where both the motion and the jump mechanism are deterministic. More recently, [92] follows [91] and proposes an estimator of $\lambda$, this time
without the assumption of a deterministic kernel. In all papers, a single trajectory of the process is fully
observed in long time, and the PDMP is supposed to be ergodic.

In [28] the authors present a method to estimate the transition kernel of a PDMP and compare it with
the Monte-Carlo estimator on a reliability problem. [4] give a recursive estimation procedure to estimate
the kernel of a PDMP from its perfect observations within a long time interval, without assumptions of
stationarity. The works presented in [7] and [9] focus on estimating PDMP-related functions rather than
inferring intrinsic process characteristics. In all those works, the statistical approach varies according to
the application domain, as each has its own specificity.

## 2.2 Impulse control for PDMPs

In this section, we delve into the control of Piecewise Deterministic Markov Processes. Broadly speaking, a stochastic control problem involves influencing the dynamics of a process through control variables or decisions to either maximize a reward or minimize a cost. The primary questions of interest include the characterization and regularity properties of the *value function* — representing the optimal reward or cost achievable across all possible decisions or decision sequences — and the existence and properties of optimal decision sequences, called *strategies*.

We begin by briefly outlining the two principal classes of control problems for PDMPs in Section 2.2.1. Our focus then shifts to impulse control problems, which we define in Section 2.2, illustrate with a detailed example in Section 2.2.3, and address in terms of resolution strategies in Section 2.2.4. Finally, we conclude this part with a discussion on the special case of control problems involving hidden information in Section 2.2.5.

### 2.2.1 Continuous versus impulse control

There are two types of control for PDMPs: *continuous control* and *impulse control*, see [45, Chapters 4 and 5]. Continuous control describes situations in which the control variable is selected at *all times* and acts on the process through the characteristics ($\Phi$, $\lambda$, and $Q$) by influencing the deterministic motion and the jump mechanism. By opposition, impulse control intervenes *punctually* on the process by moving it to a new point of the state space at some times specified by the controller.

In continuous control problems, the local characteristics of the PDMP, the flow $\Phi$, the jump intensity $\lambda$ and the transition kernel $Q$, depend on an additional control parameter that is selected by the controller at *each instant of time*. Such problems have been studied in the literature mostly under the stationarity condition by considering two strategies, the resolution of an integro-differential optimality inequality, the so called Hamilton-Jacobi-Bellman (HJB) equation (see for instance [59, 58]), or the translation into a discrete-time Markov Decision Process (MDP) on a functional space (see [1, 23, 65, 115]). Continuous control problems where the process is not required to reach stationarity have also received some attention (see e.g. [41]) but remain seldom studied.

In this paper we focus only on impulse control, where at *discrete* controller-chosen times, the process is sent to a new controller-chosen point in the state space. The simplest form of impulse control is *optimal stopping* where the controller selects a single date at which the process is stopped. Impulse control problems for PDMPs have been studied from the theoretical point of view: characterization of the value function and optimal policies through dynamic programming in e.g. [73, 38, 32] for optimal stopping and [39, 50] for the general impulse control problem. Numerical schemes to approximate the value function and compute explicit policies close to optimality have also been proposed in [39, 46].

### 2.2.2 Definition of an impulse control problem

The formal probabilistic apparatus necessary to rigorously define the impulse control problem for PDMPs is rather technical, and will not be used in the sequel, therefore, for the sake of simplicity, we only present an informal description of the problem. The interested reader is referred to [39] or [62] for a proper definition.

**Definition 3** (Impulse strategy). *A general impulse strategy S = ($\tau_n$, $\chi_n$)$_{n \ge 1}$ is a sequence of non-anticipative $\mathbb{R}_+$-valued intervention times $\tau_n$ and $\cup \subset E$-valued non-anticipative random variables $\chi_n$.*

*The trajectory of the PDMP controlled by strategy S is described in Algorithm 3. Between intervention dates, the PDMP follows its natural trajectory. At the n-th impulsion (line 12), the process is instantaneously moved to the new starting point $\chi_n$ and restarts from there.*
**Algorithm 3** Simulation of a trajectory of a PDMP controlled by strategy $S$ starting from state $x_0$ up to the $n$-th intervention

1: $x \leftarrow x_0$

2: $T \leftarrow 0$

3: **for** $n \leftarrow 1$ to $N$ **do**

4: **while** $T < \tau_n$ **do**

5: $S \sim \lambda(x)$

6: $S \leftarrow \min\{S, t^*(x)\}$

7: **if** $S < \tau_n - T$ **then**

8: $X_t \leftarrow \Phi(x, t-T)$ for $T \le t < T + S$

9: $x' \sim Q(\cdot|\Phi(x, S))$

10: $x \leftarrow x'$

11: $T \leftarrow T + S$

12: **else**

13: $X_t \leftarrow \Phi(x, t-T)$ for $T \le t < \tau_n$

14: $x' \leftarrow \chi_n$

15: $x \leftarrow x'$

16: $T \leftarrow \tau_n$

17: **end if**

18: **end while**

19: **end for**

The set of admissible strategies is denoted by $\mathcal{S}$.

**Definition 4** (Cost of an impulse strategy). The cost of strategy $S$ for a PDMP starting at $x$ at time 0 and up to infinity is defined by

$$ \mathcal{V}(S, x) = \mathbb{E}_{x}^{S} \left[ \int_{0}^{+\infty} e^{-\gamma t} c_{R}(X_{t}) dt + \sum_{n=1}^{\infty} e^{-\gamma \tau_{n}} c_{I} \left( X_{\tau_{n}}, X_{\tau_{n}^{+}} \right) \right], $$

where $\gamma \ge 0$ is a discount factor, $c_R$ is the running cost and $c_I$ is the impulse cost.

The specific assumptions that the cost functions must satisfy are omitted here, but can be found in detail in [39]. For simplicity, we present an infinite-horizon criterion. Alternative criteria are also possible, such as those with a fixed horizon $H$ or a random horizon defined by a prescribed number of jumps or impulses, as briefly discussed after Theorem 1.

**Definition 5** (Value function). The value function for a PDMP starting at $x$ at time 0 is defined by

$$ \mathcal{V}^*(x) = \inf_{S \in \mathcal{S}} \mathcal{V}(S, x). $$

In most impulse control problems for PDMPs, there is no optimal strategy as the infimum may not be reached. We define instead $\epsilon$-optimal strategies.

**Definition 6** ($\epsilon$-optimal strategy). An $\epsilon$-optimal strategy $S_\epsilon$ is an admissible strategy satisfying

$$ \mathcal{V}(S_\epsilon, x) \le \mathcal{V}^*(x) + \epsilon, $$

for all $x \in E$.

Under suitable assumptions (see again [39] for details), the value function can be characterized as the unique solution of some dynamic programming equations.

**Theorem 1** (Dynamic programming). The optimal value function $\mathcal{V}^*$ is limit of the following recursion. Set $\mathcal{V}_0^*(x) = \mathcal{V}(S_0, x)$ and $\mathcal{V}_{n+1}^*(x) = \mathcal{L}(\mathcal{V}_n^*)(x)$ for all $n \ge 0$, then for all $x$ in $E$, one has

$$ \mathcal{V}^*(x) = \lim_{n \to +\infty} \mathcal{V}_n^*(x), $$
where $\mathcal{V}(S_\emptyset, x)$ is the cost associated to the no-impulse strategy $S_\emptyset = (\tau_1 = +\infty)$ starting from $x$

$$ \mathcal{V}(S_\emptyset, x) = \mathbb{E}_x^{S_\emptyset} \left[ \int_0^\infty e^{-\gamma t} c_R(X_t) dt \right], $$

and $\mathcal{L}$ is the single-jump-or-intervention operator defined as

$$ \mathcal{L}V(x) = \inf_{t \in \mathbb{R}^+} J(MV, V)(x, t) \wedge KV(x), $$

with

$$ MV(x) = \inf_{\chi \in \mathcal{U}} \{c_I(x, \chi) + V(\chi)\}, $$

$$ J(V, W)(x, t) = \int_{0}^{t \wedge t^{*}(x)} e^{-\gamma s - \Lambda(x, s)} (c_{R}(\Phi(x, s)) + \lambda(\Phi(x, s))QW(\Phi(x, s))) ds \\ + e^{-\gamma(t \wedge t^{*}(x)) - \Lambda(x, t \wedge t^{*}(x))} V(\Phi(x, t \wedge t^{*}(x))), $$

$$ KV(x) = \int_{0}^{t^{*}(x)} e^{-\beta s - \Lambda(x, s)} (c_{R}(\Phi(x, s)) + \lambda(\Phi(x, s))QV(\Phi(x, s))) ds \\ + e^{-\gamma t^{*}(x) - \Lambda(x, t^{*}(x))} QV(\Phi(x, t^{*}(x))). $$

Operator $M$ selects the best starting point $\chi$ in the control set $\mathcal{U}$ after an intervention. Note that the control set could also depend on the position $x$ before the intervention. Operator $\inf_t J$ selects the best intervention date along the flow, while operator $K$ accounts for cases when it is best to wait for the next natural jump of the process before triggering an intervention. The iterations in the dynamic programming equation thus correspond either to natural jumps of the process or to interventions, hence the name single-jump-or-intervention operator. Dynamic programming is based on the Markov property for the canonical embedded chain, and operators can be reformulated as expectations involving the said chain, see [49].

The iterate $V_n^*$ can be interpreted as the optimal value function of the impulse control problem where at most $n$ jumps or interventions are allowed, and then strategy $S_\emptyset$ is applied. Computing the cost of the no-impulse strategy on the whole space may be demanding. Theorem 1 is also valid when replacing $\mathcal{V}(S_\emptyset, x)$ by any function $g$ such that $\mathcal{V}(S_\emptyset, x) \le g(x)$ for all $x \in E$.

Finally, we have stated the results for an infinite horizon only. If one is interested in controlling the process up to a finite time horizon $H$, one just needs to kill the process, or send it to a cemetery state, when time $H$ is reached, and define a null cost on this state to fall back to the infinite horizon framework. The discount factor $\gamma$ can be set to 0 in this context as the interval of integration is finite. Hence, the cost of strategy $S$ can typically be defined, in the fixed horizon $H$, by

$$ \mathcal{V}_{H}(S, x) = \mathbb{E}_{x}^{S} \left[ \int_{0}^{H} c_{R}(X_{t}) dt + \sum_{\tau_{n} < H} c_{I} \left( X_{\tau_{n}}, X_{\tau_{n}^{+}} \right) + c_{T}(X_{H}) \right], $$

where $c_T$ corresponds to some terminal cost function.

### 2.2.3 Running medical example as a controlled PDMP

We revisit the medical example presented in Section 2.1.4, now incorporating the additional option for the practitioner to perform surgery. This intervention entirely removes the tumor and restores the patient to a healthy state. The practitioner's objective is to prevent the patient from reaching a critical state, which occurs when the tumor size exceeds a predefined threshold, leading to death.

We add an extra mode $m = 2$ to account for the death of the patient. The state space becomes

$$ E = \bigcup_{m \in M} \{m\} \times E_{m}, $$

with $M = \{-1, 0, 1, 2\}$ and $E_{-1} = (\zeta_0, D)$, $E_0 = \{\zeta_0\} \times \mathbb{R}_+$, $E_1 = (\zeta_0, D)$ and $E_2 = \{D\}$, where $D$ is the death threshold. In this extended scenario, the biomarker dynamics are modeled by a PDMP as in Section 2.1.4 for modes $-1$, $0$ and $1$. In mode $m = 2$, the flow is constant $\Phi_2(\zeta, t) = \zeta = D$, the jump
<figure><img src="image_15.png" /><figcaption>Figure 13: **Impulse control in the medical example.** The PDMP is controlled by a threshold strategy: a surgery intervention is performed whenever the marker level exceeds the threshold $\zeta_T$. Impulsion dates are denoted by $\tau_1$ and $\tau_2$ at which the tumor is completely removed, and the marker level is reset to its nomila value $\zeta_0$.</figcaption></figure>

intensity is null $\lambda_2 = 0$ so that no more jumps are allowed, i.e. death is an absorbing state. At any time $\tau$, the practitioner may decide to intervene on the process by resetting the current state to $X_\tau = (0, \zeta_0, 0)$, representing the post-surgery state. Impulse dates can thus be chosen freely as non-anticipative dates, and the control set is a single point $U = \{(0, \zeta_0, 0)\}$.

The running cost is defined as the area under the marker curve, i.e. $c_R(\zeta) = \zeta - \zeta_0$, while the impulse cost corresponds to a fixed surgery cost, $c_I(x, x') = C$, independent of the current marker value. Intuitively, this cost structure discourages the practitioner from performing surgery prematurely when the marker level is low. Instead, the decision is delayed until the marker level reaches a threshold where the cumulative cost of waiting outweighs the surgery cost. This setup suggests a threshold-based strategy, where surgery is performed when the marker exceeds a critical level $\zeta_T$. We provide only an intuitive explanation of this result without formal proof. An example of a trajectory controlled by a strategy with intervention threshold $\zeta_T$ is shown in Figure 13.

### 2.2.4 Resolution of impulse control problems for PDMP

Optimal impulse control of PDMPs attracted attention almost as soon as the concept of PDMP was formalized by Davis. Contributions focusing only on the estimation of the value function are usually based on variational inequalities or value improvement and do not provide explicit strategies (see for instance [95, 57, 68, 56, 10]).

Explicit resolution strategies have also been proposed. In the simpler context of optimal stopping – where the decision-maker can only decide to stop the process without restarting it from a new location – the literature exhibits strategies based on either the discretization of the PDMP Markov jump kernel $Q$ [38], or on the discretization of the Markov kernel of the canonical chain [40, 51] (see Section 2.1.6 for the definition of the canonical chain). After discretization, the value function can be estimated via dynamic programming [73] using the computable discretized kernels.

These approaches have been extended for the more general case of impulse control problems. [37] considers a discretization of the state space and time to obtain uniform convergence of the discretized problem to the original one. [49] consider quantization of the kernel of the canonical chain to propose an $\epsilon$-optimal strategy. Other approaches, including iteration of a single-jump-or-intervention operator associated to the PDMP, have also successfully been proposed to construct explicit solutions in the infinite horizon context (see for instance [39, 50]).

Theoretical results on the existence of optimal strategies not relying on the discretization of kernels
have also been derived in the case of infinite horizon impulse control, but they do not provide explicit policies (see for instance [41, 42]).

All the approaches described above are numerically demanding and work for low dimensional PDMPs with limited complexity. They also require a perfect knowledge of the model and its parameters, together with a perfect observation of the process at all times (or equivalently at jumps times).

### 2.2.5 Partially observed controlled PDMPs

A very natural extension to the classical impulse control problem for PDMPs is to allow for hidden information. Suppose that the controller does not have access to the true value of the process, but only to some noisy observation of some of its components. Now decisions can only be taken based on the observations, and not on the hidden state of the process.

The literature is much scarcer for this problem. Constructive resolution approaches have been proposed in the case where jumps are still perfectly observed, but the post-jump locations are only measured through noise. In [19], the authors consider the equivalent fully observed process on the filter space and apply the quantization approach of [49] on the filtered process. [22] reduce the problem to a discrete-time Markov Decision Process (MDP) and prove the existence of optimal policies, but provide no numerical approximation of the value function or optimal strategies.

In the most challenging case where jump dates are unobserved, the authors have proposed a framework to construct an explicit strategy close to optimality in the optimal stopping context [30] and for more general impulse problems [31, 52]. Those approaches are based on a conversion of the impulse problem into a (Partially Observed) MDP, as detailed in Section 5.2.

## 3 Markov Decision Processes

Markov Decision Processes (MDPs) provide a robust framework for modeling scenarios where sequential decisions are made in the presence of uncertainty. The objective is to determine an optimal sequence of decisions that maximizes rewards or minimizes costs, while accounting for the probabilistic evolution of the system. At each step, a decision-maker observes the current state, selects an action, and receives feedback based on that action, with the system transitioning probabilistically to the next state.

MDPs belong to a class of stochastic control problems that originated in the 1950s [14, 84]. Comprehensive treatments of the subject can be found in foundational texts such as [18] and [106]. Their ability to model complex, dynamic systems has made MDPs essential across a wide range of applications, from robotics and finance to healthcare and artificial intelligence.

The versatility of MDPs in capturing stochastic dynamics and optimizing sequential decisions has led to their adoption across diverse fields, such as engineering, computer science, economics and other social sciences, making them a valuable tool for researchers and practitioners alike. In robotics, MDPs have been employed for tasks such as navigation and control [26]. Conversational agents and dialogue systems have leveraged MDPs to model and optimize human-machine interactions [96, 132]. The medical field has witnessed numerous applications of MDPs, including modeling disease progression and treatment planning for conditions like Alzheimer's disease [83], HIV [87], and diabetes [111]. MDPs have also been used in clinical decision-support systems for critical care [113].

In this section, we delve into the formal components of MDPs, including states, actions, and reward functions, while also introducing the concept of policies, which determine the action selection process (Section 3.1). Additionally, we discuss various optimality criteria in MDPs. Unless otherwise stated, we consider finite spaces for both actions and states. We then introduce several extensions of MDPs, such as the partially observed case (Section 3.2) or cases where the model is partially unknown (Section 3.3).

### 3.1 Markov Decision Process definition and control

We start with the formal definition of finite state action spaces MDPs in Section 3.1.1, then we turn to the definition of the different aspects of their optimization: policies (Section 3.1.2), evaluation criteria (Section 3.1.3) and optimal policies (Section 3.1.4). We then illustrate these concepts with a new version of our running medical example in Section 3.1.5, and briefly discuss extensions to more complex frameworks in Section 3.1.6.
<figure><img src="image_16.png" /><figcaption>Figure 14: **Graphical representation of a Markov Decision Process (MDP) with finite horizon** $H = 3$. Circles denote states, diamonds represent actions and rectangles indicate costs associated with state-action-state transitions. The dashed rectangle depicts the terminal cost.</figcaption></figure>

### 3.1.1 Generic definition

Informally, a Markov Decision Process (MDP) can be thought of as a Markov chain over a state space $S$, where transitions are influenced by actions chosen in a set $A$ by a decision-maker and evolve over an horizon of $H$ time steps. The decision process can span a finite or infinite number of time steps $H$, depending on the problem at hand.

The stochastic transitions between states are governed by a transition probability matrix $P$ (or a sequence of matrices in the case of a non-stationary MDP). The associated costs of these transitions are captured in a (possibly time-dependent) cost function $c$, which quantifies the immediate cost or reward associated with each state-action pair.

**Definition 7 (MDP).** A Markov Decision Process is defined by a tuple $\langle S, A, H, P, c \rangle$

* $S$ is the set of all possible states in the environment. In the case of discrete state spaces, $S$ may be a finite or countably infinite set.

* $A$ is the set of actions (considered finite, unless otherwise precised) that the decision-maker can take.

* $H$ is the horizon of the process, that is a positive integer (or positive infinity) counting the number of decision steps.

* $P$ describes the probabilities of transitioning, at any time step $0 \le t < H$ from a state $s_t \in S$ to another state $s' \in S$ after taking a specific action $a_t \in A$. Transition probabilities satisfy the Markov property

$$P(s'|s_0, a_0, \dots, s_t, a_t) = P(s'|s_t, a_t).$$

In the case of finite state/action spaces, transition probabilities are typically represented as transition matrices for all $a \in A$. When the MDP is non-stationary, $P = (P_t)_{t=0:H-1}$ is a family of time-indexed transition matrices.

* $c$ is the cost function that assigns a numerical cost to each state-action-state transition. $c: S \times A \times S \to \mathbb{R}$. Depending on the context, it may be worth considering a reward function instead of a cost function. As for transitions, in the case of non-stationary MDPs, $c = (c_t)_{t=0:H-1}$ is a family of cost functions.

Let us consider a finite-horizon MDP $\langle S, A, H, P, c \rangle$ and an arbitrary initial state $s_0 \in S$. At any time $0 \le t < H$, the decision-maker chooses an action $a_t \in A$ and applies it in the current state $s_t$. This moves the process to the next state $s_{t+1}$, following probability $P(\cdot|s_t, a_t)$ and yields a cost $c_t(s_t, a_t, s_{t+1})$. This process is iterated until termination when the horizon $H$ is reached. A terminal cost $C: S \to \mathbb{R}$ can be defined at the end of the process. This is illustrated on Figure 14.
**Algorithm 4** Simulation of a trajectory of a MDP controlled by policy $\pi$ starting from state $s_0$ up to the horizon $H$

1: $t \leftarrow 0$

2: $S_t \leftarrow s_0$

3: **for** $t \leftarrow 0$ to $H-1$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;$A_t \sim \pi(S_t)$

5: &nbsp;&nbsp;&nbsp;&nbsp;$S_{t+1} \sim P(\cdot|s_t, a_t)$

6: &nbsp;&nbsp;&nbsp;&nbsp;$C_t \leftarrow c(S_t, A_t, S_{t+1})$

7: &nbsp;&nbsp;&nbsp;&nbsp;$t \leftarrow t+1$

8: **end for**

### 3.1.2 MDP policies

To control an MDP, it is essential to define a *policy*, which determines the actions to be taken at each decision step based on the past *history* of the process. We now introduce the concepts necessary to define, evaluate and choose the optimal policies.

**Definition 8 (MDP history).** Let an MDP $\langle S, A, H, P, c \rangle$ be given. For any $0 < t < H$, any sequence $\langle s_0, a_0, \dots, s_{t-1}, a_{t-1}, s_t \rangle$ of states in $S$ and actions in $A$ is called (length $t$) history of the MDP. A length $0$ history is a singleton $s_0 \in S$. The set of length $t$ histories is denoted $\mathcal{H}_t$.

Then, let us define history-dependent decision rules and policies. Let $\Delta(A)$ be the set of probability distributions over $A$.

**Definition 9 (Stochastic, history-dependent decision rules).** A stochastic, history-dependent decision rule at time $t$, $\pi_t: \mathcal{H}_t \to \Delta(A)$ assigns a probability distribution over actions to every length $t$ history.

Interesting subclasses of decision rules are the deterministic and/or history-independent decision rules.

**Definition 10 (Deterministic decision rule).** A decision rule $\pi_t$ is deterministic if it assigns a unique action to every history, i.e. is a Dirac distribution over $A$.

**Definition 11 (History-independent decision rule).** A decision rule at time $t$, $\pi_t: \mathcal{H}_t \to \Delta(A)$ is history-independent if $\pi_t(h_t) = \pi_t(h'_t)$ for any pair $h_t = \langle s_0, a_0, \dots, s_{t-1}, a_{t-1}, s_t \rangle$, $h'_t = \langle s'_0, a'_0, \dots, s'_{t-1}, a'_{t-1}, s'_t \rangle \in \mathcal{H}_t$ such that $s_t = s'_t$.

Then, MDP policies are time-indexed lists of decision rules that specify the behavior of a decision-maker in an MDP.

**Definition 12 (MDP policy).** Let an MDP $\langle S, A, H, P, c \rangle$ be given. A policy for this MDP is a (potentially infinite) family of decision rules $\pi = (\pi_t)_{0 \le t < H}$. Let $\Pi_t$ denote the set of possible decision rules at time $t$ and $\Pi = (\Pi_t)_{0 \le t < H}$ denotes the full set of possible policies.

Again, specific subclasses of policies can be defined.

* $\Pi^D$ is the subset of deterministic policies, i.e. policies composed of deterministic decision rules only.

* $\Pi^M$ is the subset of Markovian policies, i.e. policies composed of history-independent decision rules only.

* $\Pi^S$ is the subset of stationary policies, i.e. Markovian policies which do not depend on $t$, these are especially useful in the case of infinite-horizon MDP.

Given a policy $\pi \in \Pi$, it is straightforward to simulate a trajectory from an MDP using Algorithm 4.

### 3.1.3 Policy Evaluation criteria

The objective in Markov Decision Processes (MDPs) is to find a policy that minimizes a cost criterion (or equivalently maximizes a reward criterion) over the entire decision horizon. The transition cost, $c(s_t, a_t, s_{t+1})$, represents the immediate penalty of taking action $a_t$ in state $s_t$ and transitioning to state $s_{t+1}$, while the terminal cost (in finite-horizon), $C(s_H)$, captures the penalty upon reaching the final state
of the process. Beyond immediate costs, we must define a criterion that captures the entire stochastic dynamics of the process, considering both the short-term effects of an action and its long-term consequences. The classical evaluation criterion is the (potentially discounted) *expected sum* of future costs.

The optimization problem associated with a cost criterion involves finding a policy $\pi$ minimizing the expected sum of costs over the set of all possible state-action trajectories, following $\pi$. We distinguish between two common cases: finite-horizon and infinite-horizon optimization criteria.

**Infinite-horizon.** In this setting, the decision-maker makes decisions indefinitely. In this case, the discounted cost criterion represents the expectation of future transition costs, discounted by a factor $\gamma \in (0, 1)$, which balances the immediate cost against future outcomes and ensures the finiteness of the sum as long as the running cost is bounded.

**Definition 13 (Discounted expected cost).** Starting from the initial state $s \in S$ and following the policy $\pi \in \Pi$ the discounted expected cost is

$$V_{\gamma}(\pi, s) = \mathbb{E}_{s}^{\pi} \left[ \sum_{t=0}^{\infty} \gamma^{t} c(S_{t}, A_{t}, S_{t+1}) \right].$$

The average cost criterion, introduced in [106], evaluates policies by comparing the limiting average (over time) of incurred expected costs.

**Definition 14 (Average expected cost).** Starting from the initial state $s \in S$ and following the policy $\pi \in \Pi$ the average expected cost is

$$V_{av}(\pi, s) = \lim_{H \to +\infty} \frac{1}{H} \mathbb{E}_{s}^{\pi} \left[ \sum_{t=0}^{H-1} c(S_t, A_t, S_{t+1}) \right].$$

**Finite-horizon.** In this setting, the decision-maker is tasked with controlling the system over a fixed number of steps, $H$. The natural objective in this case is to minimize the total expected cost:

**Definition 15 (Total expected cost).** Starting from the initial state $s \in S$ and following the policy $\pi \in \Pi$ the total cost is the sum of the costs accumulated over each step, up to the horizon $H$, plus the terminal cost associated with the final state. The expected cost can be expressed as

$$V_{H}(\pi, s) = \mathbb{E}_{s}^{\pi} \left[ \sum_{t=0}^{H-1} c_{t}(S_{t}, A_{t}, S_{t+1}) + C(S_{H}) \right].$$

Note that the above evaluation criteria apply to general history-dependent stochastic policies as well as all subclasses of policies.

### 3.1.4 Computing optimal policies

When studying an MDP, one is interested in computing an *optimal policy*, that is a policy with minimal value in every starting states.

**Definition 16 (Optimal policy).** Let MDP $\langle S, A, H, P, c \rangle$ be given. An optimal policy $\pi^*$ of the MDP, with respect to criterion $V \in \{V_{\gamma}, V_{av}, V_{H}\}$ is a policy such that:

$$V(\pi^*, s) \leq V(\pi, s), \forall \pi \in \Pi, \forall s \in S.$$

Optimal policies are not necessarily unique for a given problem. However, it is well-known (see e.g. [106]) that, in the finite state and action spaces case,

* an optimal policy exists, for any MDP and criterion $V \in \{V_{\gamma}, V_{av}, V_{H}\}$.

* an optimal stationary deterministic policy exists for $V \in \{V_{\gamma}, V_{av}\}$, in the infinite horizon case.

* an optimal (potentially non-stationary) Markovian deterministic policy exists for $V_{H}$, in the finite-horizon case.
Even though optimal policies need not be unique, they all share the same value function $V^*: S \to \mathbb{R}$, which is unique.

**Definition 17 (Value function).** Let MDP $\langle S, A, H, P, c \rangle$ be given. Let any optimal policy $\pi^*$ of the MDP, with respect to criterion $V \in \{V_\gamma, V_{av}, V_H\}$ be also given. The optimal value function (which is unique) is then defined as

$$V^*(s) = V(\pi^*, s) = \sup_{\pi \in \Pi} V(\pi, s), \forall s \in S.$$

The optimal value functions and optimal policies of an MDP with respect to the different criteria can be obtained by solving sets of non-linear equations. For the infinite-horizon discounted criterion and for the finite-horizon total cost criterion, the corresponding systems of equations are the so-called *Bellman equations* (see e.g. [14]).

**Theorem 2 (Bellman equation).**

* infinite horizon: the optimal value function $V_\gamma^*$ is the unique solution of a system of non-linear equations in variables $\{V_\gamma(s)\}_{s \in S}$, called *Bellman equations*

$$V_\gamma(s) = \min_{a \in A} \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + \gamma V_\gamma(s')), \forall s \in S.$$

Any optimal policy is such that

$$\pi^*(s) \in \arg \min_{a \in A} \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + \gamma V_\gamma^*(s')), \forall s \in S.$$

* finite horizon: the optimal value function $V^*$ is the unique solution to the following linear system of Bellman equations in variables $\{V_t(s)\}_{0 \le t \le H, s \in S}$

$$V_t(s) = \min_{a \in A} \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + V_{t+1}(s')), \forall t \in \{0, \dots, H-1\}, \forall s \in S$$

where $V_H(s) = C(s)$ defines the terminal values and $V^* = V_0$. Any optimal policy $\pi^* = (\pi_t^*)_{t=0:H-1}$ is such that

$$\pi_t^*(s) \in \arg \min_{a \in A} \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + V_{t+1}^*(s')), \forall s \in S.$$

There exists a similar defining set of equations for the average criterion, but it is omitted here due to its more technical aspects. One may refer to [106] for a definition.

We describe in Section 4 a few classical algorithms for computing optimal policies in finite and infinite horizon. Among these algorithms, classical reinforcement-learning algorithms (e.g. Q-learning introduced in [131]) make use of a related notion to the optimal-value function of an MDP: the *Q-function*.

**Definition 18 (Q-function, discounted criterion).** Let $\langle S, A, H, P, c \rangle$ be a given infinite-horizon MDP and $\pi \in \Pi$ be an arbitrary policy. The Q-function, for the infinite-horizon discounted criterion, represents the expected cost of selecting action $a \in A$ in the current state $s \in S$ and then following the policy $\pi$ forever.

$$Q_\gamma^\pi(s, a) = \mathbb{E}_s^\pi \left[ c(s, a, S_1) + \sum_{t=1}^\infty \gamma^t c(S_t, A_t, S_{t+1}) \right] \\ = \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + \gamma V_\gamma(\pi, s')).$$

Then, if $V_\gamma^*$ is the optimal value function of the MDP for the discounted criterion, $Q_\gamma^*$, the optimal Q-function of the MDP is defined as:

$$Q_\gamma^*(s, a) = \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + \gamma V_\gamma^*(s')), \forall (s, a) \in S \times A.$$

In the finite-horizon discounted case, the Q function is defined as follows.
**Definition 19 (Q-function, total cost criterion).** For any $0 \le t < H$, the Q-function at time $t$ represents the expected cost of selecting action $a \in A$ in the current state $s \in S$ and then following policy $\pi$ up to horizon $H$.

$$ 
\begin{aligned}
Q_{\{t,H\}}^{\pi}(s, a) &= \mathbb{E}_{s}^{\pi} \left[ c(s, a, S_1) + \sum_{p=t+1}^{H-1} c_p(S_p, A_p, S_{p+1}) + C(S_H) \right] \\
&= \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + V_{t+1}(\pi, s')).
\end{aligned}
$$

Then, if $\{V_{*t}\}_{t=0..H}$ form an optimal value function, the optimal Q-function is:

$$ Q_{\{t,H\}}^{*}(s, a) = \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + V_{t+1}^{*}(s')). $$

### 3.1.5 Running medical example as a finite state and action MDP

Let us consider a simplified model of a patient's disease progression using an MDP with a finite state space. It is an adaptation to the MDP framework of the example presented in Sections 2.1.2, 2.1.4 and 2.2.3.

The patient's state $s = (m, \zeta) \in S$ has two components.

* $m$ represents the patient's condition, where $m = 0$ represents remission, $m = 1$ indicates slow relapse, $m = 2$ indicates aggressive relapse, and $m = 3$ corresponds to death and is an absorbing state.

* $\zeta$ corresponds to the blood marker value, assumed to be an integer, belonging to $\llbracket 0 ; 40 \rrbracket$.

An action $a$ corresponds to the treatment choice ($a = 0$: no treatment, $a = 1$: treatment). Without treatment,

* if the patient is in remission, the marker level stays at $\zeta' = \zeta = 0$, the patient remains in a healthy condition ($m' = 0$) or transitions to a slow ($m' = 1$) or an aggressive relapse ($m' = 2$) with probabilities $p_{m'}^0$;

* during a slow relapse the biomarker $\zeta$ increases by 1 unit at each time step $\zeta' = \zeta + 1$, if the biomarker value reaches 40, the patient dies, i.e., $m' = 3$;

* during an aggressive relapse, the biomarker increases by 2 units at each time step $\zeta' = \zeta + 2$, if the biomarker value reaches 40, the patient dies, i.e., $m' = 3$.

When treatment is applied,

* if the patient is in remission, the marker level stays at $\zeta' = \zeta = 0$, and no relapse is possible;

* the biomarker decreases by 1 unit at each time step $\zeta' = \zeta - 1$, regardless of whether the relapse is slow ($m = 1$) or aggressive ($m = 2$), if the biomarker reaches 0, the patient returns to remission, i.e., $m' = 0$ and the marker value remains at $\zeta' = 0$.

At each time step, a cost is incurred based on the current state and the chosen action.

* Relapses incur a cost, with aggressive relapses ($m = 2$) having higher costs than slow relapses ($m = 1$).

* Treatment ($a = 1$) also incurs a cost.

* A terminal cost is assigned at the end of the process if the patient has died ($m = 3$).
The follow-up process is scheduled over $H = 160$ visits. The MDP associated to this example is described by

$$ 
\begin{align*}
S &= \{(0,0)\} \cup \{1,2\} \times [0; 39] \cup \{(3,40)\}, \\
A &= \{0,1\}, \\
H &= 160,
\end{align*}
$$

$$ P = \begin{cases} P(s'|(0,0),0) = p_0^0 \mathbb{1}_{\{s'=(0,0)\}} + p_1^0 \mathbb{1}_{\{s'=(1,0)\}} + p_2^0 \mathbb{1}_{\{s'=(2,0)\}} \\ P(s'|(0,0),1) = \mathbb{1}_{\{s'=(0,0)\}} \\ P(s'|(1,\zeta),0) = \mathbb{1}_{\{s'=(1,\zeta+1)\}} \mathbb{1}_{\zeta < D-1} + \mathbb{1}_{\{s'=(3,D)\}} \mathbb{1}_{\zeta = D-1} \\ P(s'|(1,\zeta),1) = \mathbb{1}_{\{s'=(1,\zeta-1)\}} \mathbb{1}_{\zeta > 1} + \mathbb{1}_{\{s'=(0,0)\}} \mathbb{1}_{\zeta = 1} \\ P(s'|(2,\zeta),0) = \mathbb{1}_{\{s'=(2,\zeta+2)\}} \mathbb{1}_{\zeta < D-2} + \mathbb{1}_{\{s'=(3,D)\}} \mathbb{1}_{\zeta \ge D-2} \\ P(s'|(2,\zeta),1) = \mathbb{1}_{\{s'=(2,\zeta-1)\}} \mathbb{1}_{\zeta > 1} + \mathbb{1}_{\{s'=(0,0)\}} \mathbb{1}_{\zeta = 1} \\ P(s'|(3,D),0) = P(s'|(3,D),1) = \mathbb{1}_{\{s'=(3,D)\}}, \end{cases} $$
$$ c(s,a) = 2a + 2 \times 1_{m=1} + 3 \times 1_{m=2}, $$
$$ C(s) = 200 \times 1_{m=3}. $$

Note that, for this example, the cost function is time-independent. The numerical resolution of this MDP is provided in the accompanying notebook with the algorithm presented in Section 4.1.1.

### 3.1.6 Generalizations

For simplicity, we defined MDPs with finite or countable state and action spaces and assumed that all actions are feasible in all states. These assumptions can be relaxed to model more complex scenarios.

**Constraints.** In many applications, specific actions may be prohibited in certain states. For example, in the medical scenario of Section 3.1.5, administering a treatment ($a = 1$) to a deceased patient ($m = 3$) should not be allowed. There are two main approaches to handle such constraints. The first one is to assign a prohibitive cost ($+\infty$ for instance) to invalid state-action pairs, preventing them from being chosen by optimization algorithms. The second one is to explicit constraints through the addition of an element $K$ to the MDP definition: $\langle S, A, H, P, c \rangle$ becomes $\langle S, A, H, K, P, c \rangle$, where $K \subset S \times A$ is called the constraints set and contains all the feasible state-action pairs. Its sections are denoted by $\{K(s), s \in S\}$, where $K(s) \subset A$ is the collection of all allowed actions in state $s$. This approach is explicit and intuitive for finite state and action spaces, requiring only that $K(s)$ be non empty for all $s \in S$. In the case of countable action spaces, it is usually required that $K(s)$ be non empty and finite for all $s$, see e.g. [80, Section 1.2] or [82, Section 2.2]. Decision rules $\pi_t$ are then restricted to select actions only from $K(s)$ for a given state or history. Such decision rules are called *admissible*, and *admissible policies* are collections of admissible decision rules.

**General state or action spaces.** MDPs can also be extended to general (uncountable) state and action spaces. While the core definitions remain similar (e.g., transition probabilities replaced by transition kernels), and in most frameworks some version of Theorem 2 still holds as well as existence of optimal policies in suitable classes of policies, these extensions introduce significant mathematical challenges.

1. In order to properly define controlled trajectories, $S_t$ and $A_t$ from Algorithm 4 must be random variables, i.e. measurable with respect to some suitably chosen $\sigma$-algebras. This first means that the state $S$ and action $A$ spaces must be endowed with $\sigma$-algebras and thus have some minimal topological structure (typically Polish spaces). This also means that *admissible* decision rules must also be measurable, and that the constraints set $K$ must contain the graph of measurable mapping from $S \times A$ onto $A$. Finally, the cost functions also need to be measurable. There exists a wide literature on easy sufficient conditions to ensure all these properties, for instance if $S$ is a subset of some $\mathbb{R}^d$ endowed with its Borel $\sigma$-algebra, and all $K(s)$ are intervals with continuous bounds in $s$, see e.g. [82, chapter 16].

2. Unbounded cost functions in general spaces may lead to undefined expectations. This issue is
resolved by assuming cost functions and transition kernels are bounded by a non-negative bounding function and using weighted norms, see [82, Chapter 16]).

3. Finally, existence of optimal policies requires some minimal regularity on the cost functions and transition kernel (typically lower or upper semi-continuity) and compactness of the set of policies (hence the use of stochastic decision rules). The interested reader is referred to [80, 23, 82] for further details.

**Semi-Markov MDPs.** As detailed in Definition 7, MDPs rely on the transition probability $P$ which verifies the Markov property: $P(s'|s_0, a_0, \dots, s_t, a_t) = P(s'|s_t, a_t)$. Similarly to the PDMP framework, it is easy to encompass semi-Markov dynamics (where the transition to a new state may depend on the time spent in the current state) by augmenting the dimension of the state space to include time in the state. In the literature, the term Semi-Markov Decision Process is also used to describe frameworks where decisions do not happen at regular intervals of time, but instead where selected actions may lead to a random amount of time before another action can be taken (this is also sometimes referred to as continuous-time MDPs). These aspects will not be discussed here, the interested reader might refer to [74] for more details, and to [90] for a formulation as an augmented MDP.

## 3.2 Partially Observed Markov Decision Processes

Partially-Observed MDPs (POMDPs) are an extension of MDPs that were first introduced by [134] within the context of control theory. In this work, we adopt a more operational perspective of POMDPs, proposed by [86], which focuses on practical algorithms for computing optimal strategies in settings with finite state, observation and action spaces. We start with a definition of POMDPs in Section 3.2.1. Then, we introduce an equivalent version as an MDP on an enlarged state space in Section 3.2.2. Finally we give an example of POMDP in our medical example context in Section 3.2.3.

### 3.2.1 POMDP definition

Here we focus on finite-horizon POMDPs.

**Definition 20 (POMDP).** A finite-horizon POMDP is a tuple $(S, A, H, P, \Omega, O, c, \omega_0)$. In addition to the components of a standard MDP introduced in Definition 7, POMDPs introduce

* a finite set $\Omega$ representing the possible observations of the system;

* an observation function

$$O: S \times A \times \Omega \to [0, 1],$$

where $O(s', a, \omega)$ gives the probability of observing $\omega \in \Omega$, after taking action $a$, resulting in successor state $s'$;

* and an initial observation $\omega_0$, describing the initial imperfect knowledge we have about the initial state.

The construction of a POMDP trajectory is illustrated on Figure 15. In a POMDP, the objective is also to compute a strategy which minimizes the expected sum of future costs. However, a key difference is that, at each time step $t$, the true state of the system $s_t$ is not directly observed. Thus POMDP histories only include actions and observations:

$$h_t = \langle a_0, \omega_1, \dots, a_{t-1}, \omega_t \rangle.$$

Similarly to solving an MDP, solving a POMDP consists of computing a *policy* assigning an action to every possible state of current knowledge of the POMDP.

**Definition 21 (POMDP history-dependent decision rule).** Consider a POMDP $(S, A, H, P, \Omega, O, c, \omega_0)$. A POMDP history-dependent decision rule $\pi_t: \mathcal{H}_t \to \Delta(A)$ assigns a probability distribution over actions in $A$ to every POMDP history $h_t \in \mathcal{H}_t$.
<figure><img src="image_17.png" /><figcaption>Figure 15: **Graphical representation of a Partially Observable Markov Decision Process (POMDP) with finite horizon $H = 3$.** Circles denote hidden states and available observations, diamonds represent actions and rectangles indicate costs associated with state-action-state transitions. The dashed rectangle depicts the terminal cost.</figcaption></figure>

**Definition 22 (POMDP policies).** Consider a POMDP $\langle S, A, H, P, \Omega, O, c, \omega_0 \rangle$. A POMDP history-dependent policy is a list of history-dependent decision rules $\pi = (\pi_t)_{0 \le t < H}$. Denote by $\Pi^o$ the set of all history-dependent policies.

Note that in practice, many algorithms seek to obtain deterministic history-dependent decision policies, i.e. lists of deterministic history-dependent decision rules $\pi_t: \mathcal{H}_t \to A$ which assign an action $a \in A$ to every POMDP history $h_t \in \mathcal{H}_t$.

A POMDP $\langle S, A, H, P, \Omega, O, c, \omega_0 \rangle$ can also be seen as an MDP $\langle \bar{S}, A, H, \bar{P}, c \rangle$ by setting

$$ \bar{S} = S \times \Omega, \\ \bar{P}(\bar{s}' = (s', \omega')|\bar{s} = (s, \omega), a) = P(s'|s, a)O(s', a, \omega'), $$

for all $\bar{s}' = (s', \omega')$ and $\bar{s} = (s, \omega)$ in $\bar{S}$. In this setting the term partially observed is more obvious. The main difference between this MPD and the POMDP setting is that in the POMPD, only policies in $\Pi^o$ are considered, and the initial value $S_0$ is randomly drawn from an initial belief $b_0$ which is a probability distribution over $S$: $b_0(s) = \mathbb{P}(S_0 = s|\omega_0)$ for all $s \in S$. This setting allows to see $\Pi^o$ as a subset of $\Pi$ for the MDP $\langle \bar{S}, A, H, \bar{P}, c \rangle$, and to directly use definitions 13 to 15 to define the cost for a POMPD policy. For any criterion $V \in \{V_\gamma, V_{av}, V_H\}$, one has

$$ V(\pi, \omega_0) = \sum_{s \in S} b_0(s)V(\pi, s). $$

**Definition 23 (Optimal policy).** Consider the POMDP $\langle S, A, H, P, \Omega, O, c, \omega_0 \rangle$. An optimal policy $\pi^*$ of the POMDP, with respect to criterion $V \in \{V_\gamma, V_{av}, V_H\}$ given in definitions 13 to 15 is a history-dependent policy $\pi^* \in \Pi$ such that:

$$ V(\pi^*, s) \le V(\pi, s), \quad \forall \pi \in \Pi^o, \forall s \in S. $$

Note that here, the expected cost $V \in \{V_\gamma, V_{av}, V_H\}$ depends on the unobserved states of the process, but must be optimized with respect only to the history.
### 3.2.2 Belief MDP

In order to evaluate a POMDP policy and compute an optimal policy, it is useful to notice that a POMDP can be modeled as another particular form of MDP over belief states. The knowledge incorporated in the history $h_t$ can be encapsulated in a belief state, $b_t$ which is a probability distribution over $S$ and represents an estimation of the (unobserved) current state. This belief state evolves from the initial belief $b_0$, based on the POMDP history $h_t \in \mathcal{H}_t$. When $h_t$ is completed with a new pair of action-observation $(a_t, \omega_{t+1})$, the belief state $b_t$ is updated as follows. Considering that in a given belief state $b_t$, an action $a_t$ is applied and leads to a new observation $\omega_{t+1}$ then

$$b_{t+1} = \tau(b_t, a_t, \omega_{t+1}).$$

where $\tau$ is a deterministic belief transition function. This function is computed through Bayesian updating:

$$ 
\begin{aligned}
b_{t+1}(s') &= \mathbb{P}(S_{t+1} = s'|b_t, a_t, \omega_{t+1}) \\
&= \frac{O(s', a_t, \omega_{t+1}) \times \mathbb{P}(S_{t+1} = s'|b_t, a_t)}{\sum_{s'' \in S} O(s'', a_t, \omega_{t+1}) \times \mathbb{P}(S_{t+1} = s''|b_t, a_t)} \\
&= \frac{O(s', a_t, \omega_{t+1}) \times \sum_{s \in S} b_t(s) P(s'|s, a_t)}{\sum_{s'' \in S} O(s', a_t, \omega_{t+1}) \times \sum_{s \in S} b_t(s) P(s''|s, a_t)} \\
&\propto O(s', a_t, \omega_{t+1}) \times \sum_{s \in S} b_t(s) P(s'|s, a_t).
\end{aligned}
$$

A way to deal with partial observability in POMDP is to consider a corresponding belief MDP, which (continuous, multi-dimensional) state space is made of the set of possible belief states of the MDP. The correspondence between POMDP and belief MDP is made explicit in the following definition.

**Definition 24 (belief MDP).** A finite-horizon POMDP $\langle S, A, H, P, \Omega, O, c, \omega_0 \rangle$ can be viewed as a belief MDP, that is an MDP $\langle \mathcal{B}, A, H, \text{Pr}, \rho \rangle$ on belief states, where:

* $\mathcal{B} \subset \Delta(S)$ is the continuous set of belief states,

* $A$ is the action space, unchanged,

* $\mathcal{P}$ is the transition function

$$ \mathcal{P}(b_{t+1}|b_t, a_t) = \sum_{\omega' \in \Omega} \mathbb{1}_{b_{t+1}=\tau(b_t, a_t, \omega')} O(s', a_t, \omega') \sum_{s \in S} b_t(s) P(s'|s, a_t). $$

* $\rho(b, a, b')$ is the cost associated to action $a$ applied in belief state $b$ resulting to $b'$

$$ \rho(b, a, b') = \sum_{s \in S} b(s) \sum_{s' \in S} b'(s') c(s, a, s'). $$

In essence, $\rho$ corresponds to the expectation of the running cost over the states following distributions given by the beliefs.

The value of a history-dependent (stochastic or deterministic) policy can now be defined as follows.

**Definition 25 (POMDP/belief MDP policy value).** Consider a finite-horizon POMDP $\langle S, A, H, P, \Omega, O, c, \omega_0 \rangle$ and the corresponding belief MDP $\langle \mathcal{B}, A, H, \mathcal{P}, \rho \rangle$, as well as a history-dependent policy $\pi$. The value (expected cost) of $\pi$ is defined as

$$ V_H(\pi, b_0) = E_{b_0}^{\pi} \left[ \sum_{t=0}^{H-1} \rho(b_t, \pi_t(h_t), b_{t+1}) \right], $$

where

$$ \rho(b_t, \pi_t(h_t), b_{t+1}) = \sum_{a \in A} \rho(b_t, a, b_{t+1}) \pi_t(h_t)(a). $$

It is well known (see e.g. [134]) that the belief MDP policy value and POMDP policy value are equal, and that optimal policies are history-dependent, deterministic in general for the POMDP version, and Markovian, deterministic for the belief MDP, which state space is finite (albeit of exponential size in the horizon) in the case of finite horizon problems [86].
### 3.2.3 Running medical example as a POMDP

Building on the MDP described in 3.1.5, we now consider that the patient's overall condition is hidden, except when the patient has died ($m = 3$), and that the biomarker is observed either through discrete or continuous noise.

**Discrete Observation Space.** Consider that the biomarker is observed through discrete noise

$$y = \zeta + \epsilon \text{ with } \epsilon \sim \mathcal{U}\{-2, -1, 0, 1, 2\},$$

where $y$ denotes the noisy biomarker. Let $z = 1_{(m=3)}$ be the death indicator. At a given time $t$, the observation of a patient's condition is $\omega = (y, z)$, and we assume that the process starts in state $s_0 = (0, 0)$. The POMDP associated with this example builds on the MDP from Section 3.1.5 and is defined by adding

$$\Omega = [-2, 42] \times \{0, 1\},$$

$$O(s', a, \omega = (y, z)) = \begin{cases} \frac{1}{5} 1_{y \in [\zeta' - 2; \zeta' + 2]; z=0} & \text{if } m' \neq 3, \\ \frac{1}{5} 1_{y \in [38; 42]; z=1} & \text{if } m' = 3, \end{cases} \\ b_0 = \delta_{(0,0)}.$$

**Continuous Observation Space.** In practice, partial observation typically involves continuous noise, and a biomarker observation is more realistically represented as

$$y = \zeta + \epsilon \text{ with } \epsilon \sim \mathcal{U}[-2, 2]$$

where $y$ still denotes the noisy biomarker. With a continuous observation space, the POMDP associated with this example builds on the MDP from Section 3.1.5 and is defined by adding

$$\Omega = [-2, 42] \times \{0, 1\},$$

$$o(s', a, \omega = (y, z)) = \begin{cases} \frac{1}{4} 1_{y \in [\zeta' - 2; \zeta' + 2]; z=0} & \text{if } m' \neq 3, \\ \frac{1}{4} 1_{y \in [38; 42]; z=1} & \text{if } m' = 3, \end{cases} \\ b_0 = \delta_{(0,0)}.$$

where $o$ is the density of the observation function $O$ with respect to the Lebesgue measure on $\mathbb{R}$.

Simulations of both the POMDP and the equivalent belief MDP for discrete observation space are provided in the accompanying notebook. An example of resolution is also presented with the python POMDP library [133] that integrates the POMDP-solvers presented in Section 4.

## 3.3 Bayesian approach to handle MDPs with unknown model

In this Section, we describe a Bayesian modeling approach of MDPs and POMDPs where the knowledge about the transition and cost functions is incomplete.

Duff et al. [61] initially proposed a Bayesian approach to represent and solve finite horizon MDPs with an unknown transition model, where transition probabilities of the MDP are modeled via a priori Dirichlet distributions. An optimal policy to solve the MDP while taking into account the uncertainty about the model is proposed. In essence, Duff shows that this optimal planning problem can be modeled as an MDP over a finite state space, where the state includes *hyperparameters* and is exponentially larger than that of the initial MDP. The implementation aspects of this framework are studied in [24] who provide empirically efficient algorithms to solve Bayes-Adaptive Markov Decision Processes (BAMDPs). An excellent review of algorithmic challenges and state of the art of the domain is given in [2]. Here we provide an overview of the BAMDP framework. We start with the definition of fully observed models in Section 3.3.1, and illustrate it with our running medical example in Section 3.3.2. Then we turn to the partially observed case in Section 3.3.3 and its medical example in Section 3.3.4.
### 3.3.1 Bayes-Adaptive MDP

In this Section we follow the didactic exposition of BAMDP and most of the notations of [110]. An archetypal BAMDP problem is when we face a (finite-horizon) MDP $\langle S, A, H, P, c \rangle$, where $S$ and $A$ are finite, the model is initially known, except for the transition probabilities $P(\cdot|s, a), \forall (s, a) \in S \times A$.

BAMDPs are embedded in a Bayesian framework, where we assume that we have a prior distribution over $P(\cdot|s, a)$. As part of the acquisition of the new state $s'$, the decision-maker updates the posterior distribution with the information provided by the transition $(s, a, s')$ when it is encountered.

Let $S = \{s^1, \dots, s^n\}$ be the finite state-space and let us define

$$p^{s,a} = (p_1^{s,a}, \dots, p_n^{s,a}),$$

where $p_i^{s,a} = P(s^i|s, a)$. The vector $p^{s,a}$ can be seen as the parameter vector of a multinomial distribution. In the Bayesian context, the conjugate of a multinomial distribution is a Dirichlet distribution, Dirichlet($\theta$), parameterized by a count vector $\theta = (\theta_1, \dots, \theta_n) \in \mathbb{N}^n$. The multinomial/Dirichlet framework makes it easy to perform Bayesian updates of the knowledge over $p$ when transitions are experienced. Indeed, assume that from a state $s$ and action $a$, we observe $K$ transitions from the true distribution $p^{s,a}$, and let $k = (k_1, \dots, k_n)$ represent the number of transitions to $s^1, \dots, s^n$ respectively (with $K = \sum_i^n k_i$). Then, a Bayesian update of a Dirichlet($\theta_1, \dots, \theta_n$) prior on $p$ leads to the posterior distribution over $p'$, equal to

$$p' = (p'_1, \dots, p'_n) \sim \text{Dirichlet}(\theta_1 + k_1, \dots, \theta_n + k_n).$$

This form of update is particularly simple. Thus, most works, starting with [61], make the assumption of independent Dirichlet priors over independent multinomial distributions $\{p^{s,a}\}_{(s,a) \in S \times A}$.

The fundamental principle of Bayes-Adaptive MDP is to incorporate model uncertainty into planning, by considering an augmented MDP, which states $\tilde{s} = (s, \theta)$ include the hyperparameters $\theta$ which can be considered as a memory of past experiences (a sufficient statistics). As in classical MDP, solving a BAMDP requires identifying the optimal decision for all possible states and hyperparameter values. Formally, a BAMDP is defined as follows.

**Definition 26 (BAMDP).** A Bayes-Adaptive MDP $\langle \tilde{S}, \tilde{A}, H, \tilde{P}, \tilde{c} \rangle$ is an MDP that can be defined from a MDP $\langle S, A, H, P, c \rangle$ with ill-known transition matrix $P$ as

* $\tilde{S} = S \times \Theta$, where $\Theta = \{(\theta_{s'}^{s,a})_{s,s' \in S, a \in A}\} \subset \mathbb{N}^{|S|^2 \times |A|}$. States $\tilde{s} = (s, \theta)$ are usually called hyper-states. Usually, an initial belief over $P$ is modeled by specifying an initial set of counts $\{\hat{\theta}_{s'}^{s,a}\}$. This initial set of counts will be assumed to verify $\sum_{s' \in S} \hat{\theta}_{s'}^{s,a} > 0, \forall s, a$, in order for the initial Dirichlet distributions to be well-defined.

* $\tilde{A} = A$. The action space is unchanged.

* $\tilde{P}$ is defined as follows: $\tilde{P}(s', \theta'|s, a, \theta)$ is the probability of transitioning from hyperstate $(s, \theta)$ to hyperstate $(s', \theta')$ under action $a$. This writes, by the chain rule:

$$\tilde{P}(s', \theta'|s, a, \theta) = \mathbb{P}(S_{n+1} = s'|S_n = s, A_n = a, \theta_n = \theta)\mathbb{P}(\theta_{n+1} = \theta'|S_n = s, A_n = a, S_{n+1} = s'\theta_n = \theta),$$
where $\mathbb{P}(S_{n+1} = s'|S_n = s, A_n = a, \theta_n = \theta) = P_\theta(s'|s, a)$ and $P_\theta$ is the transition matrix of the original MDP with parameters $p^{s,a}$ drawn from Dirichlet($\theta$), and $\mathbb{P}(\theta_{n+1} = \theta'|S_n = s, A_n = a, S_{n+1} = s'\theta_n = \theta)$ is the (deterministic) update rule described above for the hyperparameter:

$$\mathbb{P}(\theta_{n+1} = \theta'|S_n = s, A_n = a, S_{n+1} = s'\theta_n = \theta) = 1 \quad \text{iff} \quad \theta_{s'}^{s,a} = \theta_{s',a} + 1 \quad \text{and} \quad \theta_{s'}^{\hat{s},\hat{a}} = \theta_{\hat{s}',\hat{a}}, \forall (\hat{s}, \hat{a}, \hat{s}') \neq (s, a, s'). \tag{1}$$

* $\tilde{c}(s, a, s', \theta) = c(s, a, s')$, since we assume that the cost function is known.

The BAMDP $\langle \tilde{S}, \tilde{A}, H, \tilde{P}, \tilde{c} \rangle$ is a classical MDP with finite state space if $H$ is finite. In this case, the size of $\Theta$ is bounded by $H^{|S|^2|A|}$ and the size of $\tilde{S}$ by $|S| \times H^{|S|^2|A|}$. Note that, in a BAMDP, policies $\pi: \tilde{S} \to A$ choose actions as functions of the current state, $s$ and of the current knowledge about the unknown MDP transition function, summarized by the parameter vector $\theta$.

Let us make clear that the Bayes-Adaptive framework handles parameters' uncertainty in two phases.
<figure><img src="image_18.png" /><figcaption>Figure 16: Evolution of the state space of a BAMDP over a single time step. This figure illustrates a specific example with two possible states $S = \{1, 2\}$ and two possible actions $A = \{\alpha, \beta\}$. Diamonds represent actions. The parameters $\theta_{s'}^{s,a}$ correspond to the transition counts for each $(s, a, s')$ triplet. The dashed arrows indicate transition probability to hyperstates, defined according to the count matrix and the Multinomial/Dirichlet conjugacy. Here, only the case of $s_t = 1$ is shown of the Figure, but the case of $s_t = 2$ is similar.</figcaption></figure>

* The first phase is completely off-line and does not require to observe any transition of the MDP. In this phase, we use the prior knowledge about transitions and the Multinomial/Dirichlet conjugacy in order to build a (BA)-MDP with augmented state space $\tilde{S} = S \times \Theta$, as shown in Figure 16. This MDP can be solved using any MDP solution algorithm (such as those that will be mentioned in the following Sections).

* In the second phase, the decision-maker applies the optimal BAMDP policy $\tilde{\pi}: S \times \Theta \to A$ provided by the first phase to the real-life problem, in which only transitions $(s, a, s')$ drawn from the original unknown matrices $P$ are observed as shown in Figure 17. Retrieving the hyperparameter values needed to select optimal actions is straightforward since $\theta$ and $\theta'$ are transition counts and the transition $(s, a, s', \theta) \to \theta'$ is actually deterministic, as seen in Equation (1).

BAMDPs can be extended to the case where the cost function is also unknown, but we will leave this case aside. Interested readers may refer to [61] for an introduction. In [2], the authors provide an overview of the literature on algorithmic approaches to solve BAMDP, as well as a complexity analysis of the full problem. Among the algorithms, we will briefly mention in Section 4.3.3 a Reinforcement Learning approach to solve BAMDP, based on Monte-Carlo planning.

### 3.3.2 Running medical example as a BAMDP

Going back to the MDP medical example presented in Section 3.1.5, we now suppose that the parameters $p_0^0, p_1^0$ and $p_2^0$ are unknown, i.e. that the probability transitions from a health condition state (0,0) to any of the states (0,0), (1,0) or (2,0) in the absence of treatment are unknown. Recall that all other transitions are deterministic.

To model the uncertainty on these transitions corresponding to $P(\cdot|s = (0,0), a = 0)$, we use the standard Multinomial/Dirichlet approach

$$ (p_0^0, p_1^0, p_2^0) \sim \text{Dirichlet}(\theta_0^0, \theta_1^0, \theta_2^0). $$

With each new observation, we update the a posteriori distribution of the parameters of the transition matrix $\theta' = \theta + \Delta$, where $\Delta = e_{m'}$ is the vector in $\mathbb{N}^3$ that has a value of 1 at the position corresponding to the observed new mode $m'$, while other entries are 0. Let us illustrate the update process with an example. Suppose that at time $t$, the patient is in state $s = (0,0)$ and no treatment ($a = 0$) is applied. At this stage, the hyperparameter vector $\theta$ is, for instance, $\theta = (5, 1, 0)$. At the next step, the patient transitions into a slow relapse state. This observation corresponds to the state-action-state sequence
<figure><img src="image_19.png" /><figcaption>Figure 17: **Graphical representation of BAMDP with finite horizon H = 3.** Circles denote states and hyperparameters, dashed rounded rectangles denote hyperstates, diamonds represent actions and rectangles indicate costs associated with state-action-state transitions. The dashed rectangle depicts the terminal cost.</figcaption></figure>

$s = (0, 0), a = 0, s' = (1, 0)$. In this case, the parameter $\theta_1^0$ of the Dirichlet distribution is updated, and $\Delta$ is equal to $e_2 = (0, 1, 0)$. The hyperparameter vector $\theta$ is updated as follows

$$ \theta' = (5 \ 1 \ 0) + (0 \ 1 \ 0) = (5 \ 2 \ 0). $$

The BAMDP associated with this example builds on the MDP from Section 3.1.5. It is defined by adding

$$ \begin{aligned} \tilde{S} &= S \times [0; H]^3, \\ \tilde{A} &= A, \\ \tilde{c} &= c, \end{aligned} $$

$$ \tilde{P}(s' = (m', \zeta'), \theta' | s, a, \theta) = \begin{cases} P(s'|s, a) \mathbb{1}_{\theta' = \theta}, & \text{if } (s, a) \neq (0, 0, 0), \\ \frac{\theta_{m'}^0}{\theta_0^0 + \theta_1^0 + \theta_2^0} \mathbb{1}_{\theta' = \theta + e_{m'}}, & \text{if } (s, a) = (0, 0, 0). \end{cases} $$

Resolution of the BAMDP is provided in the accompanying notebook with algorithms presented in Section 4.3.3. Results are illustrated through trajectory simulations.

### 3.3.3 Bayes-Adaptive POMDP

In [110], the author proposed an extension of BAMDP to the partially-observed case. In this *Bayes-Adaptive POMDP* model (BAPOMDP), the state of the process is hidden and the only available information is given through observations $\omega \in \Omega$. In addition, the observation model is parameterized, similarly to the transition model, and the values of the parameters are updated through interactions with the system.

**Definition 27 (BAPOMDP).** A BAPOMDP is an extension of a BAMDP given by $\langle \tilde{S}, A, H, \tilde{P}, \Omega, O, c, \bar{b}_0 \rangle$, where

$$ \bullet \ \bar{S} = \tilde{S} \times \Psi = S \times \Theta \times \Psi, \text{ where} $$

$$ \Psi = \left\{ (\psi_{s', \omega}^a)_{a \in A, s' \in S, \omega \in \Omega} \right\} \subset \mathbb{N}^{|S| \times |A| \times |\Omega|} $$

is an observation count function. The coordinate $\psi_{s', \omega}^a$ counts the number of times $\omega$ has been obtained when action $a$ has led to $s'$. Note that it is not precisely observed since $s'$ is not observed.
• $\tilde{P}$ is defined using Bayes rule and independence properties as $\tilde{P}(s', \theta', \psi', \omega|s, \theta, \psi, a)$ is defined as (by Bayes rule and independence properties):

$$ 
\begin{aligned}
\tilde{P}(s', \theta', \psi', \omega|s, \theta, \psi, a) \\
&= \mathbb{P}(S_{t+1} = s'|S_t = s, A_t = a, \theta_t = \theta) \times \mathbb{P}(\theta_{t+1} = \theta'|S_t = s, A_t = a, S_{t+1} = s', \theta_t = \theta) \\
& \quad \times \mathbb{P}(\psi_{t+1} = \psi'|\psi_t = \psi, A_t = a, S_{t+1} = s', \omega_{t+1} = \omega) \times \mathbb{P}(\omega_{t+1} = \omega|A_t = a, S_{t+1} = s', \psi_t = \psi), \\
&= \tilde{P}(s', \theta'|s, a, \theta) \times \mathbb{P}(\psi_{t+1} = \psi'|\psi_t = \psi, A_t = a, S_{t+1} = s', \omega_{t+1} = \omega) \times \mathbb{P}(\omega_{t+1} = \omega|A_t = a, S_{t+1} = s', \psi_t = \psi).
\end{aligned}
$$

As for the Dirichlet state transition probability, the observation probability is defined as

$$ \mathbb{P}(\omega_{t+1} = \omega|A_t = a, S_{t+1} = s', \psi_t = \psi) = \frac{\psi_{s', \omega}^a}{\sum_{\omega' \in \Omega} \psi_{s', \omega'}^a}, $$

and the transition to $\psi'$ is also deterministic

$$ \mathbb{P}(\psi_{t+1} = \psi'|\psi_t = \psi, A_t = a, S_{t+1} = s', \omega_{t+1} = \omega) = 1 \quad \text{iff} \quad \psi_{s', \omega}^{a} = \psi_{s', \omega}^{a} + 1 \quad \text{and} \quad \psi_{\hat{s}, \hat{\omega}}^{a} = \psi_{\hat{s}, \hat{\omega}}^{a}, \forall (\hat{a}, \hat{s}, \hat{\omega}) \neq (a, s', \omega). $$

• Finally, the initial belief state is defined over $\bar{S}$ as: $\bar{b}_0(s, \theta_0, \psi_0) = b_0(s)$ and $\bar{b}_0(s, \theta, \psi) = 0$ if $(\theta, \psi) \neq (\theta_0, \psi_0)$. The initial counts $(\theta_0, \psi_0)$ are assumed to be known.

As for BAMDPs, the BAPOMDP model $\langle \bar{S}, A, H, \tilde{P}, \Omega, O, c, \bar{b}_0 \rangle$ is a new POMDP over an extended finite state space (provided that $H$ is finite). In [110] the authors provide exact (for very small problems) and approximate solution algorithms, for solving BAPOMDP, see Section 4.

### 3.3.4 Running medical example as a BAPOMDP

Going back to the BAMDP example presented in section 3.3.2, we now consider a noisy measurement $y$ of the biomarker as introduced in section 3.2.3: $y = \zeta + \epsilon$ with a multinomial noise $\epsilon \in \{-2, -1, 0, 1, 2\}$ with unknown parameters $p_i = p_{s'=(m', \zeta')}, \omega=y \mathbb{1}_{\zeta'-y=i}$ as $\mathbb{P}(\epsilon = i)$ independent from the action. We assume that the $p = (p_{-2}, p_{-1}, p_0, p_1, p_2)$ follows a Dirichlet distribution. At the initial stage, the parameters of the observation model are given by

$$ p \sim \text{Dirichlet}(\psi_2, \psi_{-1}, \psi_0, \psi_1, \psi_2), $$

with $(\psi_i)_{-2:2} \in \Psi \subset \mathbb{N}^5$.

The update of the parameters for $\psi$ is similar to that of the parameters for $\theta$ presented in section 3.3.2. With each new observation $(s', a, \omega)$, we update the posterior distribution of the parameters of the transition matrix, so that $\psi' = \psi + \Delta$ is a counting matrix, where $\Delta = e_i$ is the vector in $\mathbb{N}^5$ that has a value of 1 at the position corresponding to $i = \zeta' - y$ and 0 elsewhere, for $s' = (m', \zeta')$ and $\omega = y$.

This BAPOMDP is therefore equivalent to a POMDP on a larger state-space. As in the classic POMDP formalism (see Section 3.2.1), this uncertainty can be represented using regular POMDP formalism presented in Section 3.2.1, by maintaining (and updating) a belief over these count vectors.

## 4 Algorithms for computing exact or approximate optimal policies in MDPs

The exact computation of optimal policies is often challenging due to combinatorial explosion, the curse of dimensionality and/or the continuous nature of many real-world problems. To address these challenges, various methods have been developed to approximate optimal policies. This section explores these methods, beginning with exact solution techniques and transitioning to advanced approximation methods, such as simulation based methods or Reinforcement Learning (RL).

Exact solution techniques involve targeting the optimal policy with convergence guarantees, often through smart exploration of the whole set of admissible policies. While these methods are theoretically appealing, they quickly become computationally infeasible for large-scale problems due to the curse of dimensionality. To overcome these limitations, approximation techniques provide a practical alternative. Many of these approaches rely on simulations to strike a balance between precision and computational efficiency. Monte Carlo Planning (MCP) is a notable example, leveraging random sampling to evaluate
<caption>**Algorithm 5 Backward Iteration Algorithm**</caption>

1: **for** $s \in S$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;$V_H(s) = C(s)$
3: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
4: **for** $t$ from $H-1$ to $0$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**for** $s \in S$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$V_t(s) \leftarrow \min_{a \in A} \sum_{s' \in S} P(s'|s, a)(c_t(s, a, s') + V_{t+1}(s'))$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\pi_t^*(s) \leftarrow \arg \min_{a \in A} \sum_{s' \in S} P(s'|s, a)(c_t(s, a, s') + V_{t+1}(s'))$
8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
9: **end for**

and refine policies. Techniques like Monte Carlo Tree Search (MCTS) excel in sequential decision-making problems by simulating state-action trajectories and updating estimates based on averaged outcomes. Beyond Monte Carlo methods, approaches such as approximate dynamic programming and policy iteration with sampling have been proposed. These methods offer additional practical solutions for complex control problems but fall outside the scope of this paper. The interested reader is referred to [27] for further details on this topic.

Reinforcement Learning provides a powerful framework for approximating value functions and deriving near-optimal policies, especially in cases where exact methods are impractical. RL methods focus on learning optimal policies through iterative interactions with the environment, making them particularly effective in scenarios with partially unknown or highly complex models.

This section offers a high-level overview of these approaches, covering their theoretical underpinnings and practical applications without delving into algorithmic details. Specifically, Section 4.1 introduces dynamic programming-based algorithms, which provide exact solutions; Section 4.2 discusses Monte Carlo Planning algorithms, which are designed to handle problems of higher complexity, and finally Section 4.3 presents RL approaches, tailored for frameworks with partially unknown models.

## 4.1 Dynamic programming algorithms

MDPs with finite state and action spaces and explicit models may be solved exactly via dynamic programming. This iterative method is based on Bellman equations, stated in Theorem 2. It calculates the optimal value function and policy by defining and solving simpler embedded subproblems. This section extensively uses notations and definitions from Section 3.1 and describes dynamic programming for finite-horizon and then infinite horizon problems. We start with backward iteration in Section 4.1.1, then introduce value iteration and policy iteration in Sections 4.1.2 and 4.1.3 and conclude with a discussion on limitations and extensions in Section 4.1.4.

### 4.1.1 Backward iteration

The backwards induction algorithm given in Algorithm 5 consists of systematically solving decision-making problems over a fixed number of time steps by working backwards from the final stage. Unlike the infinite-horizon case, where the goal is to optimize decisions over an indefinite period, in the finite-horizon case, decisions are optimized over a fixed number of stages. The process begins with the definition of the terminal value function, where no further decisions can be made. It thus corresponds to the terminal cost function. For each preceding time step, the value function is updated using Bellman recursion from Theorem 2, which calculates the expected cumulative cost based on the current state, available actions, immediate costs, and the expected value of future states. This backward induction process continues until we reach the initial time step $t=0$, at which point the value function for each state has been determined for every stage. Once the value function has been calculated for all stages, a Markovian optimal policy can be derived by selecting, at each time step, the action that minimizes the expected cost for each state. This approach ensures that the resulting policy is optimal for the finite-horizon problem, as it takes into account both immediate and future costs at each decision point. In finite-horizon MDPs, the optimal decision rule generally varies with the time step $t$, making the optimal policy non-stationary. The accompanying jupyter notebook provides an implementation of the backward iteration algorithm for the MDP version of our running medical example given in Section 3.1.5.
## Algorithm 6 Value Iteration Algorithm

**Require:** Discount factor $\gamma \in (0, 1)$, convergence tolerance $\epsilon$

1: **for** $s \in S$ **do**

2: $V_0(s) \leftarrow 0$

3: **end for**

4: $n \leftarrow 0$

5: $\delta \leftarrow +\infty$

6: **while** $\delta > \epsilon(1 - \gamma)\gamma$ **do**

7: **for** $s \in S$ **do**

8: $V_{n+1}(s) \leftarrow \min_{a \in A} \sum_{s' \in S} P(s'|s, a)(c(s, a, s') + \gamma V_n(s'))$

9: **end for**

10: $\delta \leftarrow \max_{s \in S} |V_{n+1}(s) - V_n(s)|$

11: $n \leftarrow n + 1$

12: **end while**

13: **for** $s \in S$ **do**

14: $\pi^*(s) \leftarrow \arg \min_{a \in A} \sum_{s' \in S} P(s'|s, a)(c(s, a, s') + \gamma V_{n+1}(s'))$

15: **end for**

### 4.1.2 Value iteration

The value iteration algorithm given in Algorithm 6 is a dynamic programming method used to solve MDPs in infinite-horizon problems. The algorithm updates the values $V(s)$ at each iteration, for every states $s \in S$, converging to the optimal value function $V^*$. It begins with arbitrary initial values $V_0(s), \forall s \in S$ and improves them by repeatedly applying the Bellman optimality equation given in Theorem 2. The stopping criterion is based on the maximum change in the value function. Once the algorithm stops it returns a value function $V_n$, from which an optimal policy $\pi^*$ can be derived (even though $V_n$ may not equal $V^*$). Note that the optimal policy obtained this way is stationary. The convergence of the algorithm to the optimal policy is guaranteed for discounted factors $0 < \gamma < 1$ and sufficiently small $\epsilon$, see [14].

### 4.1.3 Policy iteration

The policy iteration algorithm given in Algorithm 7 is closely related to the value iteration algorithm described in Section 4.1.2. While value iteration focuses on iterating over value functions, policy iteration directly updates the policies themselves. After initialization with an arbitrary initial policy, the algorithm operates by iteratively improving the policy as follows.

1. *Policy Evaluation:* this step computes the value of each state under the current policy by solving the associated system of linear equations or iteratively estimating the value function.

2. *Policy Improvement:* this step updates the policy by deriving a *greedy* policy with respect to the current value function. This involves selecting actions that maximize the expected return for each state.

These steps are repeated, with the newly updated policy being evaluated and improved in turn. The process continues until the policy converges to the optimal policy, which remains unchanged in subsequent iterations. Policy iteration is guaranteed to compute the optimal policy in a finite number of steps, see [84]. The number of iterations required is polynomial in the number of states and actions of the MDP, making it an efficient algorithm for solving MDPs.

### 4.1.4 Limitations and extensions

The three algorithms discussed above rely on dynamic programming equations, which represent one of the most efficient methods for exploring the policy space. While dynamic programming is guaranteed to yield an optimal policy in theory, it often becomes computationally infeasible for real-world applications. This limitation arises from the curse of dimensionality, a phenomenon identified by [14], where the computational complexity grows exponentially with the dimensionality of the state and action spaces, especially if those are continuous.
**Algorithm 7 Policy Iteration Algorithm**

**Require:** Discount factor $\gamma \in (0, 1)$, initial policy $\pi_0$

1: **for** $s \in S$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp; $\pi(s) \leftarrow \pi_0(s)$

3: **end for**

4: **repeat**

5: &nbsp;&nbsp;&nbsp;&nbsp; $V \leftarrow$ solution of the following linear program in variables $(V(s), s \in S)$ \{Policy Evaluation\}

$$V(s) = \sum_{s' \in S} P(s'|s, \pi(s)) (c(s, \pi(s), s') + \gamma V(s')), \quad \forall s \in S$$

6: **Update greedy policy $\pi$** \{Policy Improvement\}

$$\pi(s) \leftarrow \arg \min_{a \in A} \sum_{s' \in S} P(s'|s, a) (c(s, a, s') + \gamma V(s')), \quad \forall s \in S$$

7: **until** $\pi$ does not change

One approach to mitigate this complexity is to employ piecewise discrete policies, enabling the application of dynamic programming techniques (see, for example, [16]). Discretization simplifies the problem by drastically reducing the dimensionality of the state space, making computations more manageable. However, this reduction comes at the cost of precision, often resulting in suboptimal policies (see [53]). Additionally, the process of discretizing a state space is itself challenging. It can be time-consuming, difficult to implement effectively, and prone to introducing inaccuracies, especially as the dimensionality of the problem increases.

Given these limitations, alternative approaches such as approximation methods and RL algorithms – explored in the following sections – are often preferred. These methods provide practical solutions for handling high-dimensional discrete or continuous spaces more effectively, frequently outperforming discretized dynamic programming in terms of both policy quality and computational efficiency.

For simplicity, the remainder of this discussion will focus on finite-horizon problems to streamline the analysis and facilitate the exploration of these alternative methods.

## 4.2 Monte-Carlo simulations-based algorithms

For finite-horizon MDPs with a known starting state, solving the MDP involves exploring a *min-expectation* tree. The root of this tree corresponds to the initial state, $s_0$. Each action $a_0 \in A$ generates a successor action node. Each action node $a_0$ further branches into successor state nodes for all $s_1$ such that $P(s_1|s_0, a_0) > 0$. This process continues across time steps, with the tree’s leaves representing the final states $s_H$.

Solving the MDP using exact dynamic programming involves fully constructing this tree and applying backward induction. Backward induction labels each state node $s_t$ with its optimal value $V_t^*(s_t)$, which represents the return obtained by following the optimal policy $\pi^*$ from $s_t$ until the horizon $H$.

However, backward induction is computationally intensive because it requires exploring the entire tree. Even storing the complete tree is infeasible in many cases, though depth-first exploration with pruning of dominated actions can help mitigate this storage issue.

To address these computational challenges, Monte Carlo Planning methods, and in particular the family of Monte Carlo Tree Search (MCTS) algorithms, offer an efficient alternative. These heuristic-free approaches, such as the Upper Confidence bound applied to Trees (UCT) algorithm [89], focus on exploring the tree through simulations. MCTS algorithms prioritize exploration of the most promising parts of the tree—areas that show high potential or require further evaluation to refine their estimated values. By doing so, they avoid exhaustive exploration, making them well-suited for large or complex MDPs (see, e.g., [43]).

There exist many variants of MCTS [21, 135], however most of them iteratively update a search tree (which is not exactly the same as the dynamic programming tree) until a computation budget is exhausted, at which point an action, computed for the root node, is returned. A tuple $(s, a, \nu, \rho)$ is attached to each node $v$. Component $v.s = s \in S$ is a state, $v.a = a \in A$ is the action that led to state
<figure><img src="image_20.png" /><figcaption>Figure 18: **The four steps of Monte-Carlo Tree Search.** A child node is selected (here node 5), then the tree is expanded with a new child node (here node 9) and a trajectory is simulated from this new node with default policy. Finally, the tree is backpropagated to update the Q-values.</figcaption></figure>

$v.s$, $v.\nu = \nu \in \mathbb{N}$ is the number of visits to node $v$ following the current edge (i.e. the number of times $v.a$ was applied in $v'$.s and led to $v.s$, if $v' = Par(v)$ is the parent node of $v$), $v.\rho = \rho \in \mathbb{R}$ is an estimate of the value function $V^*(v.s)$ at state $v.s$. Initially, the tree consists of a single root node $v_0$ with attached tuple $(v_0.s = s_0, v_0.a = \emptyset, v_0.\nu = 1, v_0.\rho = 0)$. At any step, the algorithm loops over the following four successive substeps to update the current tree, see Figure 18 for an illustration.

1. **Selection.** Starting from the root node, child nodes are selected recursively.

* If $v$ is the current node, we check that it has a child node $v^a$ for every $a \in A$, i.e. if all available actions have been tried in state $v.s$. Else, Selection phase is over (condition I).

* Then, an action $a$ is selected, using a *tree policy* and a successor state $s' \sim P(\cdot|v.s, a)$ is simulated. If there exist $v'$ such that $v'.s = s'$ and $v'.a = a$, then $v'$ is selected and the Selection phase goes on from $v'$. Else, Selection phase is over (condition II).

2. **Expansion.** When Selection phase is over, the tree is expanded. A new child node of $v$ (selected node at the end of the Selection phase) is added. It is obtained

* either (if condition I was met) by playing a yet unplayed action $a \in A$ and simulating $s' \sim P(\cdot|v.s, a)$,

* or (if condition II was met) corresponding to the newly simulated state $s'$.

Let $v'$ be this node and $(v'.s = s', v'.a = a, v'.\nu = 1, v'.\rho = 0)$ the attached tuple.

3. **Simulation.** When a node $v'$ has been added in the expansion phase, a single trajectory is simulated, starting from $v'.s$, using a *default policy* (for example, successive actions may be selected uniformly at random). Let $C$ be the total (discounted in case of a discounted MDP) cost collected along the simulated trajectory. Then, $v'.\rho \leftarrow C$.

4. **Backpropagation.** The Q-values (see definitions 18 and 19) of the ancestor nodes of $v'$ are updated recursively: For every ancestor nodes $v$ of $v'$, selected in backward order,

$$v.\rho \leftarrow \min_{a \in A} \frac{\sum_{\tilde{v} \in Succ(v), \tilde{v}.a=a} \tilde{v}.\nu \times (c(v.s, a, \tilde{v}.s) + \gamma \tilde{v}.\rho)}{\sum_{\tilde{v} \in Succ(v), \tilde{v}.a=a} \tilde{v}.\nu}.$$

Once the Backpropagation step is over, the algorithm goes back to the Selection step, until the computation budget is over. At this point, an action $a^*$ is returned for the initial state $s_0$:

$$a^* \leftarrow \arg \min_{a \in A} \sum_{v' \in Succ(v_0), v'.a=a} \frac{\sum_{\tilde{v}' \in Succ(v), \tilde{v}'.a=a} \tilde{v}'.\nu \times (c(v.s, a, \tilde{v}'.s) + \gamma \tilde{v}'.\rho)}{\sum_{\tilde{v}' \in Succ(v), \tilde{v}'.a=a} \tilde{v}'.\nu}.$$
<figure><img src="image_21.png" /><figcaption>Figure 19: **Reinforcement Learning schematics.** The learning process involves a decision-maker observing a state, taking an action that affects the state, and receiving feedback through costs.</figcaption></figure>

An MCTS algorithm instance is determined by the choice of a *tree policy* and a *default policy*. Typically, the tree policy is derived from a UCT criterion, balancing *exploitation* of the best action so far and *exploration* of under-used actions. The default policy is often uniformly random. However many choices of tree policies and default policies are possible, leading to different algorithms. The reader may refer to [135] for a recent overview.

Note that, contrarily to exact dynamic programming, Monte-Carlo planning does not require complete knowledge of the MDP model, (P, c), but only the faculty to simulate $P(\cdot|s, a)$ and to record experiences (s, a, s', c). Even though the backpropagation step uses costs $c(v.s, a, v.s')$, these costs have been encountered during the expansion step and may be added to the tree data structure. This feature is common with RL approaches, described in the next section.

## 4.3 Reinforcement-learning algorithms

Reinforcement Learning (RL) is a computational framework that enables decision-makers to learn optimal mappings from situations to actions in order to minimize a numerical cost. In RL, the environment represents everything beyond the decision-maker's direct control. In the context of MDPs, the environment corresponds to the mathematical structure that defines state transitions and costs, as described in Section 3.1. Unlike traditional learning methods, RL does not require complete knowledge of this structure. Instead, it allows decision-makers to learn directly through interactions with the environment.

The RL process involves a cycle of interaction: the decision-maker observes the current state, selects an action, and receives feedback in the form of a cost, which reflects the quality of its decision (see Figure 19). This action also influences the next state of the environment. Through repeated iterations, the decision-maker refines its behavior, ultimately developing a policy that optimizes long-term costs, as outlined in Section 3.1.

One of the fundamental challenges in RL is balancing the trade-off between exploration and exploitation. To minimize costs, an RL decision-maker must exploit actions that have proven effective in the past while also exploring new actions to uncover potentially better strategies. Striking this balance is essential to prevent the decision-maker from becoming stuck in suboptimal behaviors and to enable continuous policy refinement through exploration.

The process of learning optimal actions in RL is approached from multiple perspectives, each characterized by distinct methodologies and challenges. This section begins with *model-free methods*, which do not rely on a model of the environment's dynamics. Within this category,

* *Value-based* approaches focus on estimating the value function to derive optimal policies (Section 4.3.1);

* *Policy-based* methods directly search for the optimal policy within the policy space (Section 4.3.2).

Following model-free approaches, we delve into *model-based methods* in Section 4.3.3, which involve constructing an explicit model of the environment to predict future states and associated costs.

By structuring the discussion in this manner, we aim to provide a comprehensive overview of the key strategies used in RL, emphasizing the strengths and limitations of each approach. This section draws
**Algorithm 8 Q-learning**

**Require:** discount factor $\gamma > 0$, tolerance $\epsilon > 0$, learning rate $\alpha \in (0, 1)$, arbitrary initial values for $Q(s, a), \forall s \in S, a \in A$

1: **for** each episode **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;Initialize $s$ randomly

3: &nbsp;&nbsp;&nbsp;&nbsp;**for** each step of episode **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Choose $a \leftarrow \arg\max_{a' \in A} Q(s, a)$ with probability $\epsilon$ or uniformly in $A$ with probability $(1 - \epsilon)$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;($\epsilon$-greedy choice)

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Apply action $a$, observe $c, s'$

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$Q(s, a) \leftarrow Q(s, a) + \alpha [c + \gamma \min_{a' \in A} Q(s', a') - Q(s, a)]$

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s \leftarrow s'$

8: &nbsp;&nbsp;&nbsp;&nbsp;**end for**

9: **end for**

heavily on the seminal work by [125], which offers an extensive survey of RL methodologies and their applications.

Recent advancements in deep learning (e.g. [93, 70]) have significantly enhanced RL capabilities, particularly in addressing challenges posed by high-dimensional state and action spaces. These developments have given rise to deep Reinforcement Learning (deep RL) algorithms, which are discussed throughout this section. These algorithms show considerable promise for advancing the field.

For practitioners, practical implementation can begin by encoding their environment using the (python) *Gymnasium* library introduced by [20]. This library integrates seamlessly with the library *RLlib* [97], which offers a suite of pre-implemented deep RL algorithms ready for use.

### 4.3.1 Model-free approaches: value-based algorithms

Model-free value-based RL is a key area within artificial intelligence that focuses on learning optimal policies by estimating value functions, without the need for an explicit model of the environment. This approach is especially valuable in complex environments where creating an accurate model is either difficult or impractical.

**Tabular algorithms.** We begin with the simplest cases, where the state and action spaces are small enough for the value functions to be approximated using arrays. These methods are known as tabular algorithms, see [125].

One of the foundational algorithms in this category, *Q-learning* was introduced by [130]. It is an *off-policy* method that allows a decision-maker to learn a policy different from the one it follows during training. This distinction enables the decision-maker to learn the optimal policy while adhering to a more exploratory behavior policy. In *Q-learning*, the decision-maker interacts with the environment and updates its *Q*-values based on received costs and the minimum expected sum of future costs, as defined by the Bellman equation. The iterative process is continued until a stopping criterion is reached. [131] show that provided the action state is finite, the *Q*-values converge to the optimal action values almost surely. Through this process, the decision-maker identifies the best possible action for any given state, as outlined in Algorithm 8.

*Q-learning* is widely recognized for its simplicity and effectiveness in the case of finite state and action spaces. However, [131] also show that it can be computationally expensive and slow to converge, particularly in large state spaces.

In [112], the authors introduce the State-Action-Reward-State-Action (SARSA) algorithm, which is also model-free and value-based. The name SARSA reflects the sequence of events the algorithm uses to update its knowledge. Unlike *Q-learning*, SARSA updates its *Q*-values based on the action actually taken by the decision-maker, rather than the optimal action, making it an *on-policy* algorithm in contrast to *Q-learning*'s off-policy approach. [119] show that the SARSA algorithm is guaranteed to converge if every state-action pair is visited infinitely often, but its convergence can be slow, especially in environments with large state spaces. Another alternative to *Q-learning* is Temporal-Difference (TD) learning, introduced by [124], which combines Monte Carlo methods with dynamic programming.

Over the last few decades, extensions to these algorithms have been developed to address their limitations. For instance, an *on-policy* *Q-learning* extension was created to achieve better computational
<figure><img src="image_22.png" /><figcaption>Figure 20: **Deep Q-Network algorithms.** Deep Q-Network (DQN) algorithms send the current state of the environment to a neural network, which calculates the *Q* values for each possible action in the MDP. The algorithm then selects the action that minimizes these *Q*-values. Once the action has been chosen, the environment moves to the next state and the associated cost is calculated. This cost is used to update the loss function, which measures the difference between the predicted *Q*-values and the actual values. Finally, the network weights are adjusted by backpropagation to improve future predictions.</figcaption></figure>

complexity bounds [121], perform in multi-agent environments [116], and address the issue of overestimated action values [77].

**Function approximation techniques.** While tabular algorithms and their extensions work well in finite and relatively small state spaces, many real-world problems involve larger or even continuous state spaces. In cases where a tabular representation of the *Q*-function is impractical, [126] proposed to approximate the *Q*-function using function approximation techniques. Here, the value function is not represented as a table but as a parameterized function with a weight vector $\mathbf{w} \in \mathbb{R}^d$. The approximate value of a state $s$ given weight vector $\mathbf{w}$ is denoted by $\hat{v}(s, \mathbf{w})$. The approximation $\hat{v}$ can take various forms such as a linear function, a decision tree, or others. [125] extended the theoretical results for methods using function approximation to cases involving partial observability.

Deep RL extends function approximation techniques by using a multi-layer artificial neural network build the approximate the value function $\hat{v}$. This approach provides a powerful method for function approximation, albeit without theoretical guarantees [125]. Many Deep RL algorithms are built upon foundational tabular methods. A notable example is the Deep Q-Network (DQN), introduced by [100] and illustrated in Figure 20. Subsequent advancements, such as *Double DQN* [78], *Dueling DQN* [129], and *Rainbow* [81], have further enhanced the efficiency and performance of DQN-based approaches.

As depicted in Figure 20, when dealing with partially observable environments, the neural network receives observations rather than full state information. In such cases, the optimal policy depends on the entire history of observations, as discussed in Section 3.2.2. Addressing this challenge, [99] introduced *Recurrent Replay Distributed DQN* (R2D2), which uses recurrent neural networks (RNNs) instead of convolutional networks. RNNs are capable of retaining information from past observations, making them particularly well-suited for partially observable Markov Decision Processes (POMDPs).

These algorithms collectively demonstrate the flexibility and potential of model-free value-based RL methods in tackling complex decision-making challenges across various domains.

### 4.3.2 Model-free approaches: policy-based algorithms

While value-based RL algorithms generally learn a *Q*-function, from which a (nearly-) optimal policy can be obtained, policy-based algorithms directly learn a policy from experiments. The emblematic dynamic programming policy-based algorithm is *policy iteration*, which alternates phases of *evaluation* of the current (deterministic) policy and *improvement* of the current policy (see Section 4.1).

**Policy gradient approaches.** When the model of the MDP is unknown (RL context) and/or when state and action spaces are too large or even continuous, policy iteration is no longer applicable. In this context, *Policy Gradient approaches* (see e.g. [54, 117]) suggest to use continuously parameterized policies
and replace policy improvement steps with *gradient descent* steps, to improve the current policy (defined by its continuous parameters). With parameterized policies, the value of the policy in an arbitrary step can be considered itself as a function of the parameters of the policy. As such, the gradient of the value function can be estimated using stochastic policy gradient algorithms (the value function in state *s* is an expectation of the sum of costs obtained along possible trajectories).

More precisely, model-free *policy gradient* algorithms exploit the *gradient theorem* from [123] which states that

$$ \nabla_{\theta} V(\pi_{\theta}, s_0) = \mathbb{E} \left[ \nabla_{\theta} \log \pi_{\theta}(A|S) \tilde{Q}(S, A) \right], $$

where the expectation is taken with respect to $S \sim \mu_{s_0}^{\pi_{\theta}}$ and $A \sim \pi_{\theta}$, $\mu_{s_0}^{\pi_{\theta}}$ is the discounted occupation measure of the Markov chain followed by starting in $s_0$ and applying policy $\pi_{\theta}$, $\pi_{\theta}(\cdot|\cdot)$ is any differentiable parameterized representation of a stochastic policy (with parameter vector $\theta$) and $\tilde{Q}$ is an estimation of the *true* $Q$-function (which can be computed, e.g. by simulation). This policy gradient expression makes it easy to estimate the gradient of the current policy by simulation, since the gradient of $\pi_{\theta}$ may be known explicitly, and simulations can be used to estimate $\tilde{Q}$ as well as the global expectation. Note also that policy-gradient algorithms naturally handle continuous state and action spaces.

One limitation of these approaches is the need to arbitrarily choose a parameterized policy, something that may not always be ideal. Moreover, estimating the gradient through simulations often requires substantial computational resources, even if the model is known. The challenge becomes more significant when the model is unknown, as the approach then depends on datasets of observed MDP trajectories.

**Parametrized policy approaches.** On the other hand, policy iteration approaches are highly versatile and capable of addressing remarkably complex MDP problems. This is especially true with recent advancements handling neural network representations of policies. Proximal Policy Optimization (PPO, [114]) is a popular example of a deep RL algorithm with policy iteration. In PPO, instead of approximating the value function as in deep $Q$-networks, the policy itself is parameterized by a neural network. This network directly associates states (or observations in the case of partial observability) with actions, optimizing the policy to minimize expected costs. In PPO, the algorithm alternates between collecting data from the environment using the current policy and updating the policy on the basis of the data collected. As with DQN, PPO can be extended to handle partial observability by incorporating recurrent neural networks instead of feedforward networks. This allows the policy to retain past observations and improve decision-making in environments with hidden states, such as POMDPs.

Unlike DQN, which generates a deterministic policy by selecting the action with the highest $Q$ value, PPO generates a stochastic policy. A stochastic policy is one in which the decision-maker selects actions according to a probability distribution over the available actions, rather than always choosing a fixed action for a given state. This stochastic nature allows the decision-maker to explore more diverse behaviors, which makes it particularly useful in environments where there is uncertainty or partial observability, as it can prevent the policy from getting stuck in sub-optimal solutions by occasionally trying different actions.

### 4.3.3 Model-based algorithms

This class of algorithms operates in an iterative loop: the decision-maker begins by collecting data to update the current model of the environment's dynamics. Using the refined model, the decision-maker learns and optimizes its policy within a simulated environment. Once a new policy is derived, the decision-maker applies it in the real world, gathering additional data. This new data is then used to further refine the model. The cycle of data collection, model updating, and policy optimization repeats continuously, enabling the decision-maker to improve the performance over time. This iterative process is depicted in Figure 21.

Atkeson et.al [3] show that this provides better results than model-free approaches with less data. For the sake of simplicity, we only focus on Bayesian Reinforcement Learning (BRL) in this survey. However, readers interested in model-based methods are encouraged to explore the comprehensive survey on model-based RL provided in [98], as well as the state-of-the-art DreamerV3 algorithm proposed by [76], a model-based deep RL algorithm.

Introduced in [101], BRL is a model-based RL approach that leverages methods from Bayesian inference to model uncertainty in the environment's dynamics. The Bayesian framework, described in
<figure><img src="image_23.png" /><figcaption>Figure 21: **Model-based Reinforcement Learning.** Data is collected to update the environment’s model knowledge, which is then used to learn and optimize the policy in a simulated environment. This updated policy is applied, yielding new data.</figcaption></figure>

section 3.3.1, allows for principled incorporation of prior knowledge and can help in balancing the exploration-exploitation tradeoff effectively.

A Bayesian approach to dynamic programming was introduced by [122], where a model is sampled from the posterior distribution and dynamic programming is applied to the sampled model. To improve accuracy, the model is periodically resampled to account for new information. While effective, this method can be computationally slow, especially in complex environments. In [61], a method is proposed that utilizes Finite State Controllers (FSCs) to represent Bayesian policies in a more compact form. FSCs are graph-based structures where the nodes represent memory states and the edges correspond to transitions between states, defined by the tuple (s, a, s', c). The algorithm then searches the space of possible FSCs to identify the optimal controller, effectively reducing the complexity of finding the best policy. The Bayesian Exploration Exploitation Tradeoff in LEarning (BEETLE) algorithm [104] builds upon the PERSEUS algorithm from [120], originally developed for solving continuous POMDPs. The key innovation of BEETLE is to reformulate the BAMDP as a POMDP that operates over the space of states and transition functions. This approach leverages the strengths of POMDP techniques to address the complexities inherent in BAMDPs.

Due to the high complexity associated with model-based BRL, most approaches have been limited to very small state spaces, typically comprising only 10 to 20 states. As the number of states increases, a large amount of data is required to accurately learn a robust model. In addition, planning methods become difficult to implement in larger state spaces, as they require consideration of the full range of possible posterior distributions, which significantly complicates the decision-making process.

Online tree search methods have become a promising solution for solving BAMDPs by evaluating potential future action sequences from the current state, building a decision tree to guide actions. Unlike offline methods, online approaches focus on finding effective policies for the current belief state, making them more feasible for complex problems. A notable example is Bayesian Sparse Sampling, introduced in [128], which non-uniformly explores the tree of future states and actions, reducing computational complexity by concentrating on the most promising paths. Another well-known approach is MCTS algorithms [21] (see Section 4.2), which balance exploration and exploitation through random sampling.

Building on MCTS, [72] developed Bayesian Adaptive Monte Carlo Planning (BAMCP) to extend tree search for BAMDPs, incorporating Bayesian updates to inform action exploration. Despite improvements over simpler methods, BAMCP remains computationally demanding, especially in large state and action spaces, and lacks function approximation, limiting its use in very large or continuous state spaces.

In BRL, significant progress has been made with algorithms designed to handle continuous state spaces and complex environments. BayesianExplore (BE) [88] is a standout algorithm, leveraging posterior sampling of value functions for efficient exploration. Another key advancement is Deep Bayesian Reinforcement Learning (Deep BRL), which has been particularly effective in dealing with continuous state spaces. Methods such as BRL-Deep Sparse Sampling (DSS) from [71] combine deep learning with sparse sampling to approximate value functions in high-dimensional continuous spaces where other methods often falter.
# 5 Controlled PDMP as member of the MDP family

In this section, we demonstrate how impulse control problems for PDMPs can be reformulated within the MDP framework. We also outline how this shift in perspective opens up new avenues for addressing problems that have so far remained unresolved within the PDMP framework.

For example, as discussed in Section 2.2.2, formulating a rigorous approach to impulse control problems for PDMPs involving hidden jumps and noisy observations presents significant challenges, particularly in defining admissible strategies (see, for instance, [1] or [42, Section 1.1]). The MDP framework offers a robust structure for overcoming these challenges and for properly formalizing such problems.

Here, we show how a fully observed impulse control problem for PDMP can be formulated as an MDP (Section 5.1), and then turn to the partially observed version (Section 5.2).

## 5.1 Controlled PDMP as MDP

We start by encompassing the notion of impulse for PDMP in the notion of action for MDPs in Section 5.1.1, then formally state a sub class of impulse control problems for PDMPs as MDPs in Section 5.1.2, and finish in Section 5.1.3 with a new version of our running medical example fitting this extended framework.

### 5.1.1 Impulse control strategies and MDP policies

As detailed in Sections 2.2.2 and 3.1, the notions of impulses for PDMP and actions for MDP slightly differ. On the one hand, in impulse control, the decision-maker typically has two tasks

1. choose the next impulse date,

2. at an impulse date, select the point from which to restart the process.

On the other hand, in the MDP framework, the decision-maker takes decisions at each stage, and decisions act on the dynamics of the process by changing the probability of observing new states through the controlled transition kernel $P$.

To unify both frameworks, we choose to restrict the choice for the impulse strategies of the PDMP decision-maker. Other less restrictive approaches involving more general and complex state spaces are possible, see e.g. [62] and will be briefly discussed at the end of this section. Our main aim here is on the one hand to keep the exposition as simple as possible, and on the other hand to use finite action spaces MDPs in view of their numerical resolution. We also only consider a finite horizon in the sequel, for the same reasons.

The first restriction concerns the control set. In the impulse control framework for PDMPs, the decision-maker chooses the new starting point of the process after an impulse, which means they may choose to change either the mode or the Euclidean variable of the process or both. In the sequel, we only allow the decision-maker to change the mode of the process. In this view, we introduce a mode-augmented PDMP framework in which we consider two types of modes, one corresponding to the decision-maker-chosen modes, the other to the modes the decision-maker cannot act on. Formally, we decompose the discrete state space $M$ of a PDMP into a finite product set $L \times M$, where modes in $L$ can be decision-maker chosen while modes in $M$ cannot.

**Definition 28** (Mode-augmented controlled PDMP). A mode-augmented controlled Piecewise Deterministic Markov Process is a PDMP defined by the tuple $(E, \Phi, \lambda, Q)$, where

* The state space $E$ has the specific form

$$E = \{(\ell, m, x), \ell \in L, m \in M, x \in E_m^\ell\},$$

for some finite product mode set $M = L \times M$, and where $E_m^\ell = E_{m=(\ell,m)}$ is some Borel subset of $\mathbb{R}^{d_m}$.

* The flow $\Phi$ just satisfies the conditions from definition 1, and we write

$$\Phi(x, t) = (\ell, m, \Phi_m^\ell(x)),$$

for all $x = (\ell, m, x) \in E$.
• The jump intensity just satisfies the conditions from definition 1 and we write

$$ \lambda_m(x) = \lambda_m^\ell(x), $$

$$ \text{for all } x = (\mathbf{m} = (\ell, m), \mathbf{x}) \in E. $$

• The jump kernel $Q$ satisfies the conditions from definition 1 with the additional constraint that $Q$ cannot change the value of $\ell$, as $\ell$ is intended to be chosen by the decision-maker. We also set

$$ Q(\cdot|x) = Q_m^\ell(\cdot|x), $$

$$ \text{for all } x = (\ell, m, \mathbf{x}) \in E. $$

From this rigorously defined mode-augmented PDMP, we can further drop the regime $\ell \in L$ from the state $x \in E$ and include it in the action instead to better fit the standard MDP notation. Denote

$$ E_m = \bigcup_{\ell \in L} E_m^\ell, \quad E_M = \{(m, \mathbf{x}), m \in M, \mathbf{x} \in E_m\}. $$

We also adapt notations for the kernels $P_{t_{n+1}-t_n}$ associated with the skeleton chain defined in Section 2.1.6 as follows

$$ P_{t_{n+1}-t_n}^{\ell}(B|m, \mathbf{x}) = P_{t_{n+1}-t_n}(\{\ell\} \times B|\ell, m, \mathbf{x}) = \mathbb{P}(X_{t_{n+1}} \in \{\ell\} \times B|X_{t_n} = (\ell, m, \mathbf{x})), \quad (2) $$

$$ \text{for any Borel set } B \subset E_M \text{ and } x = (\ell, m, \mathbf{x}) \in E. $$

The second restriction concerns the impulse dates. In the MDP framework, decision dates are often deterministic and predefined. In contrast, within the impulse control framework for PDMPs, determining the timing of the next impulse is typically part of the decision-making process and is allowed to depend on the current state of the process. This introduces a key distinction: for a fixed real-time duration $H$ of the impulse control task, different strategies may result in varying numbers of *impulse dates*, meaning different numbers of decisions.

To reconcile this difference, and to keep finite action spaces, we impose a maximum number $H$ of impulses, for example by restricting impulse dates to a finite time grid. For simplicity, we assume that decision dates are limited to a subset $\mathbb{T}$ of a uniform time grid $0, \delta, \dots, H\delta$ for some $\delta > 0$. This assumption can be relaxed if needed. Keeping track of both the control horizon $H$ and the maximum number of decisions $H$ induces additional technicalities. We thus introduce an artificial decision, denoted $\tilde{d}$, which is triggered when the real-time horizon $H$ is reached. This decision will be applied repeatedly as many times as needed to align with the decision horizon $H$, transitioning the process to an absorbing cemetery state $\Delta$. This approach necessitates the introduction of constraints on the set of actions, as decision-makers are not permitted to control the process beyond the time horizon. To formalize this, we incorporate constraints into the decision framework using constraint sets $K$ as mentioned in Section 3.1.6.

In summary, an action $a$ for the MDP corresponding to the controlled PDMP will consist in both a decision-maker-chosen mode $\ell \in L$ and a time until the next decision $r \in \mathbb{T}$, to encompass both aspects of impulses for PDMPs. As mentioned above, other formalisms closer to the original impulse problem for PDMPs are possible. They involve more complex MDPs, with possibly continuous action spaces. In our approach, the strongest restriction is that the delay between consecutive impulses belongs to a fixed grid $\mathbb{T}$, and thus is not allowed to depend on the current state of the process, whereas in the original PDMP framework, decision dates are stopping-times, which allows much more flexibility. This flexibility may be retrieved by defining larger action and constraints spaces. The MDP framework can also increase the flexibility of the impulse control problem for PDMPs as MDP actions are allowed to have a random effect, whereas a PDMP impulse fixes the new starting point of the process.

### 5.1.2 Formalism of the equivalent MDP

We can now formally consider impulse control for a PDMP in an MDP framework. In order to avoid confusion between the continuous time index $t \in \mathbb{R}_+$ of the PDMP and the discrete time index $t \in \mathbb{N}$ of the MDP stages, we rename the latter with indices $n \in \mathbb{N}$.

**Definition 29** (Controlled PDMP as an MDP). Let $\langle E, \Phi, \lambda, Q \rangle$ be a mode-augmented controlled PDMP. The corresponding MDP is defined by the MDP $\langle S, A, H, K, P, c, C \rangle$ with the following characteristics.
* The state space is

$$S = E_M \times \{0, \delta, \dots, (H-1)\delta\} \cup \{\Delta\},$$

where $\Delta$ is a cemetery state, and a time variable $t_n \in \{0, \delta, \dots, H\delta\}$ is added to keep track of the continuous-time horizon.

* The horizon $H$ is finite and equals $H/\delta$ (which is supposed to be an integer).

* The action space is

$$A = (L \times \mathbb{T}) \cup \{\tilde{d}\},$$

where $\tilde{d}$ is an artificial action to deal with the continuous-time horizon of the PDMP.

* The constraints set $K \subset S \times A$, is such that its sections $K(s) = \{a \in A; (s, a) \in K\}$ satisfy

$$K(\Delta) = \{\tilde{d}\}, \\ K(s = (x, t_n)) = L \times (\mathbb{T} \cap \{\delta, \dots, H\delta - t_n\}).$$

This ensures that, in the cemetery state, the only action allowed is $\tilde{d}$, and that no decision can send the process beyond the time horizon $H = H\delta$. We assume that $\delta \in \mathbb{T}$ so that the constraints sets are never empty.

* The transition kernel $P$ is the transition kernel $P_{t_{n+1}-t_n}$ associated with the skeleton chain of the mode-augmented PDMP, where the grid $(t_0, \dots, t_H)$ corresponds to the decision dates, combined with the deterministic increase of time. More precisely, if $s = (m_{t_n}, x_{t_n}, t_n)$, and $a = (\ell, r)$, then $t_{n+1} = t_n + r$ and for any Borel subset $B$ of $E_M$, one has

$$P(s' = (x' \in B, t_{n+1}) | s = (x, t_n), a = (\ell, r)) = P_r^{\ell}(B|x) \mathbb{1}_{t_{n+1}=t_n+r},$$

where $P^\ell$ is defined in Eq. (2). Additionally, the artificial action $\tilde{d}$ sends the process to the cemetery state

$$P(\{\Delta\} | s, \tilde{d}) = 1,$$

for all $s \in E$.

* The cost function $c$ should in principle be defined by

$$c(s, a, s') = c(s = (x, t_n), a = (\ell, r)) = \mathbb{E} \left[ \int_0^r c_R(X_t) dt \mid X_0 = (\ell, x) \right] + \tilde{c}_I(s, a),$$

where $\tilde{c}_I(s, a) = c_I(X_{t_n} = (\ell_{t_n}, x), X_{t_n^+} = (\ell, x))$ if $a = (\ell, r)$, $s = (x, t_n)$, and the decision is taken to move $X_{t_n} = (\ell_{t_n}, x)$ to $X_{t_n^+} = (\ell, x)$. In practice, such a cost function is generally intractable, and alternative cost functions close to $c$, depending only on the values of $s, a, s'$ are considered instead. We provide an example in Section 5.1.3. The terminal cost $C$ is equal to $c_T$.

One of the main advantages of the MDP framework, is that it is now very easy to rigorously define decision rules $\pi_n$ as measurable functions from $S$ onto $A$ such that $\pi_n(s) \in K_n(s)$ for all $s$ and $n$, and to define admissible policies, as exposed in Section 3.1.2. The controlled trajectory of the MDP following policy $\pi = (\pi_n)_{0 \le n < H} \in \Pi$ is defined recursively in Algorithm 9. Line 4 corresponds to Algorithm 2 for the mode-augmented PDMP starting from state $x = (\ell_n, X_n)$ at time 0 and up to time horizon $r_n$. Note that the cemetery state $\Delta$ ensures that all trajectories have the same length $H$, even if they do not have the same number of actual decisions ($a \ne \tilde{d}$).

Then the total expected cost of policy $\pi \in \Pi$ starting at $s_0 \in S$ is

$$V_H(\pi, s_0) = \mathbb{E}_{s_0}^{\pi} \left[ \sum_{n=0}^{H-1} c(S_n, A_n, S_{n+1}) + C(S_H) \right],$$
**Algorithm 9** Simulation of a trajectory of a controlled PDMP seen as an MDP controlled by policy $\pi = (\pi_t)_{0 \le t < H}$ starting from state $s_0$ up to the horizon $H = \delta H$

1: $S_0 \leftarrow s_0 = (x_0, 0)$

2: **for** $n$ from 0 to $H-1$ **do**

3: $A_n = (\ell_n, r_n) = \pi_n(S_n)$

4: $S_{n+1} \sim P(\cdot|S_n = (X_n, T_n), A_n) = (P_{r_n}^{\ell_n}(\cdot|X_n), T_n + r_n)$

5: **end for**

defined in Section 3.1 is equivalent to $V_H(S, x)$ defined in Section 2.2.2, and the impulse control problem corresponds to the optimization problem

$$V(s_0) = \inf_{\pi \in \Pi} V_H(\pi, s_0).$$

To obtain a formal equality $V(x_0) = V(s_0 = (x_0, 0))$ is possible but requires a perfect match between the cost functions and the sets of admissible policies. In particular, admissible policies for the controlled PDMP must satisfy the same restrictions as those introduced in Section 5.1.1. This will not be further discussed here as we chose the MDP framework specifically to avoid delving into the formal definition of admissible policies or cost functions for PDMPs.

### 5.1.3 Running medical example as a controlled PDMP formalized as an MDP

We now consider a more interesting version of our running medical example as a controlled mode-augmented PDMP and then specify the corresponding MDP.

**Controlled PDMP** We consider again a model for the dynamics of the marker of a patient as a PDMP as described in Section 2.1.4, with very different decision from those proposed in Section 2.2.3. We now consider that impulses correspond to the choice of the next visit date to the medical center and of the treatment to be applied until this next visit. To do so, we consider a new mode and time augmented PDMP. Its state space is defined from the modes sets $L = \{0, 1\}$ and $M = \{0, 1, 2, 3\}$ and the state spaces

$$E_0^\ell = \{\zeta_0\} \times [0, H]^2, \quad E_1^\ell = E_2^\ell = [\zeta_0, D) \times [0, H], \quad E_3^\ell = \{D\} \times [0, H],$$

for all $\ell \in L$. In this setting,

* $\ell \in L$ corresponds to the treatment applied, where $\ell = 0$ indicates no treatment, and $\ell = 1$ represents the standard treatment.

* $m \in M$ still corresponds to the overall condition of the patient ($m = 0$: sound, $m = 1$: main disease, $m = 2$: therapeutic escape, $m = 3$: death),

* $\zeta \in [\zeta_0, D]$ still corresponds to the marker value, $\zeta_0$ being the nominal value for a sound patient, and we further assume that if the marker reaches level $D$, the patient dies,

* $u \in [0, H]$ (in mode $m = 0$) corresponds to the time since the last change of overall condition $m$.

Note that we deal with both a mode-augmented and time-augmented (in modes $m = (0, 0)$ and $m = (1, 0)$) PDMP.

Decision dates correspond to visits to the medical center, during which the marker is measured, and a new treatment is selected and applied until the next visit. At each visit, the practitioner must decide on two aspects: the interval until the next visit, $r \in \mathbb{T} = \{15, 30, 60\}$ days ($\delta = 15$), and the treatment to be administered $\ell \in L$. Thus, at each decision point, the practitioner selects a decision $a = (\ell, r)$. The patient's follow-up period is denoted by $H$, marking the real-time horizon of the decision process. The local characteristics of the PDMP are defined as follows, depending on the treatment $\ell$.

* If no treatment ($\ell = 0$) is applied, in the disease conditions ($m = 1$ or $m = 2$), the marker level grows exponentially and can reach the death level $D$ (boundary jump to $m = 3$) in finite time. The death level is absorbing. In the therapeutic escape condition $m = 2$, no other change of state is possible. In the main disease condition $m = 1$, a random jump to therapeutic escape is also possible. In the remission condition $m = 0$, the marker level remains at its nominal value $\zeta_0$ and the patient may randomly jump
<figure><img src="image_24.png" /><figcaption>Figure 22: Medical example as a mode-augmented PDMP. Example of a trajectory from a mode-augmented PDMP, with a strategy including only visits every 15 days. Background colors indicate patient condition mode, while colors of the trajectory indicate treatment modes (black: $\ell = 0$, green: $\ell = 1$). At time $t = 0$ the process starts with action $a = (0, 15)$; at time $\tau_1$ the decision-maker starts treatment, i.e. $a = (1, 15)$, at time $\tau_2$ the decision-maker stops the treatment, i.e. $a = (0, 15)$.</figcaption></figure>

to any of the two disease conditions $m \in \{1, 2\}$. Formally, the characteristics of the PDMP for $\ell = 0$ are thus

$$ 
\left\{ 
\begin{aligned}
\Phi_0^0(\zeta_0, t) &= \zeta = \zeta_0, \\
\Phi_1^0(\zeta, t) &= \zeta e^{v_1^0 t}, \\
\Phi_2^0(\zeta, t) &= \zeta e^{v_2^0 t}, \\
\Phi_3^0(\zeta, t) &= \zeta = D,
\end{aligned} 
\right. 
\quad 
\left\{ 
\begin{aligned}
\lambda_0^0(\zeta, u) &= f_{\to 1}^0(u) + f_{\to 2}^0(u), \\
\lambda_1^0(\zeta) &= f_{\to 2}^0(\zeta), \\
\lambda_2^0(\zeta) &= 0, \\
\lambda_3^0(\zeta) &= 0,
\end{aligned} 
\right. 
\quad 
\left\{ 
\begin{aligned}
Q_0^0(B|\zeta_0, u) &= \frac{f_{\to m'}^0(u)}{f_{\to 1}^0(u) + f_{\to 2}^0(u)} 1_B(m', \zeta_0), \\
Q_1^0(B|\zeta) &= 1_{\{\zeta < D\}} 1_B(2, \zeta) + 1_{\zeta = D} 1_B(3, D), \\
Q_2^0(B|D) &= 1_B(3, D),
\end{aligned} 
\right. 
$$

for any Borel subset $B$ of $E_M$, where $f_{\to 1}^0$ and $f_{\to 2}^0$ are intensity functions corresponding respectively to a jump from mode 0 to mode 1 and 2. For instance, they can be chosen as Weibull distributions as in Section 2.1.4.

• Under treatment ($\ell = 1$), in remission ($m = 0$) the marker level remains at its nominal value and the patient may only randomly jump to the therapeutic escape condition $m = 2$. In the disease condition $m = 1$, the marker level decreases exponentially and can reach the nominal level $\zeta_0$ ($m = 0$) in finite time or randomly jump to the therapeutic escape condition $m = 2$. In the therapeutic escape condition $m = 2$, the marker level grows exponentially and can reach the death threshold $D$ ($m = 3$) in finite time, no other change of state is possible. Again, the death state is absorbing. Formally, the characteristics of the PDMP for $\ell = 1$ are thus

$$ 
\left\{ 
\begin{aligned}
\Phi_0^1(\zeta_0, t) &= \zeta = \zeta_0, \\
\Phi_1^1(\zeta, t) &= \zeta e^{-v_1^1 t}, \\
\Phi_2^1(\zeta, t) &= \zeta e^{v_2^1 t}, \\
\Phi_3^1(\zeta, t) &= \zeta = D,
\end{aligned} 
\right. 
\quad 
\left\{ 
\begin{aligned}
\lambda_0^1 &= f_{\to 2}^1(u), \\
\lambda_1^1 &= f_{\to 2}^1(\zeta), \\
\lambda_2^1 &= 0, \\
\lambda_3^1 &= 0,
\end{aligned} 
\right. 
\quad 
\left\{ 
\begin{aligned}
Q_0^1(B|\zeta_0, u) &= 1_B(2, \zeta_0), \\
Q_1^1(B|\zeta = \zeta_0) &= 1_B(0, \zeta_0, 0), \\
Q_1^1(B|\zeta \neq \zeta_0) &= 1_B(2, \zeta), \\
Q_2^1(B|D) &= 1_B(3, D),
\end{aligned} 
\right. 
$$

for any Borel subset $B$ of $E_M$, where again $f_{\to 2}^1$ is the intensity function corresponding to a jump from mode 0 to mode 2. Figure 22 illustrates a trajectory from this PDMP with a given strategy including treatment and next decisions dates.
**Formulation as an MDP** We assume that visits take place at discrete dates $t_n$, where the time lapse between two visits can be 15, 30 or 60 days. At most $H = \frac{H}{15}$ visits can occur. Moreover, decision-related constraints appear. The last visit must take place at the end $H$ of the follow-up (We assume $H$ is a multiple of 15). The associated MDP is $\langle S, A, H, K, P, c, C \rangle$. Its characteristics are defined as follows.

* The state space $S$ corresponds to the state of a patient $x \in E_M$, a counter of time since the beginning $t$ and the cemetery state $\Delta: S = E_M \times \{0, \delta, \dots, (H-1)\delta\} \cup \{\Delta\}$,

* The action space $A$ is expressed either as a decision pair $a = (\ell, r)$, where $r \in \mathbb{T} = \{15, 30, 60\}$ is the delay until the next visit and $\ell \in L = \{0, 1\}$ is the treatment choice or as the artificial action $a = \tilde{d}: A = L \times \mathbb{T} \cup \{\tilde{d}\}$.

* The horizon $H$ is equal to $H/15$ decision dates.

* The constraints $K(s) \subset A$ is the space of admissible decisions in state $s$,

$$ \mathbb{K}(s) = \begin{cases} \tilde{d} & \text{if } s = (m, x, t) \text{ and } m = 3 \text{ or } s = \Delta, \\ L \times (\mathbb{T} \cap \{\delta, \dots, H\delta - t\}) & \text{otherwise.} \end{cases} $$

* The transition function $P$ can be written as a combination of the PDMP flow, jump intensity and Markov kernel. Detailed analytic formulas are omitted here, but notice that the MDP transition function is set according to the PDMP parameters.

* The cost function $c$ as introduced in Section 5.1.2 is approximated by

$$ c(s, a, s') = (\zeta - \zeta') \times r + C_V + C_\ell, $$

where $a = (r, \ell)$, $\zeta$ is the second coordinate of $s$ and $\zeta'$ that of $s'$ (if both $s$ and $s'$ differ from $\Delta$), i.e. the running cost is approximated by the histogram estimator of the area under the curve, and the impulse cost corresponds to a visit cost $C_V$ and a cost of treatment $C_\ell$. The terminal cost is $C_D > 0$ if the patient has died, 0 otherwise.

Figure 23 illustrates a trajectory of the PDMP controlled as an MDP, with perfect observations on visit dates.

## 5.2 Partially observed controlled PDMP as POMDP

While properly defining an impulse control problem for a hidden PDMP is very challenging, the framework of POMDPs allows to deal with partially hidden MDPs. Transforming an impulse control problem for partially observed PDMPs into a POMDP is done in a similar manner as in the perfect observation case, and again greatly simplifies the problem statement.

Here we assume that observations are only available at the decision dates $t_n$, and are given by $Y_n = F(X_{t_n}) + \epsilon_n$ with $(\epsilon_n)$ real-valued independent and identically distributed random variables with density $f$ independent from the PDMP and the actions, and $F$ a real-valued link function (to simplify the exposition. Both the noise and link function could be multi-dimensional). Only partial information about the process $X$ is obtained through $Y$, and we hence consider $X_{t_n}$ as hidden. We also assume that the time since the beginning of the follow-up is observed through a variable $w_n = t_n$. The objective is now to control the PDMP through policies that may only rely on the observations available up to each decision date.

We first state this problem as a POMDP in Section 5.2.1, then extend our running medical example to this more realistic context in Section 5.2.2, and finally discuss the computation of belief states originating from PDMPs in Section 5.2.3.

### 5.2.1 Formulation as a POMDP

We keep the restrictions on the impulses from Section 5.1.1, and now give a POMDP formulation of the partially observed impulse control problem for a PDMP.

**Definition 30** (Partially observed controlled PDMP as a POMDP). Let $\langle E, \Phi, \lambda, Q \rangle$ be a mode-augmented controlled PDMP under partial observations. The corresponding POMDP is defined by $\langle S, A, H, P, \Omega, O, K, c, C, b_0 \rangle$ with the following characteristics.
<figure><img src="image_25.png" /><figcaption>Figure 23: **Trajectory of a PDMP controlled as an MDP.** Dots indicate the perfect observations obtained only at visit dates. Here again, background colors indicate the patient condition mode, while colors of the dots indicate the treatment mode. In this example, visits are made every 15 days (i.e. $r = 15$) and the treatment changes twice.</figcaption></figure>

* The state space $S$, the action space $A$, the transition kernel $P$, the cost functions $c$ and $C$ and the horizon $H$ are the same as those of Definition 29.

* The observation space is $\Omega = (\mathbb{R} \times [0, H]) \cup \{\Delta\}$. We will denote $\omega = (y, t) \in \Omega$.

* The observation function $O$ corresponds to the observation function of the PDMP, which may or may not depend on the decision $a$. Here, it has density

$$o(s' = (x', t'), a, \omega = (y, t)) = f(y - F(x'))1_{t=t'},$$

with respect to the Lebesgue measure on $\mathbb{R}$.

* The constraints set $K \subset \Omega \times A$, satisfy the same conditions as in Definition 29 and may only depend on the observations.

* The initial belief $b_0$ is the filter at time $t_0$ associated to the partially observed PDMP, and will be discussed more thoroughly in Section 5.2.3.

Here again, the impulse control for the partially observed PDMP is equivalent to its POMDP counterpart, where both the state space $S$ and the observation space $\Omega$ are continuous. It is then straightforward to consider the associated belief MDP, which can take a particular meaning in the PDMP context, as explained in Section 5.2.3.

### 5.2.2 Running medical example as a partially observed controlled PDMP formalized as a POMDP

We consider the same example as in Section 5.1.3, but now assume that the patient condition $m$ and time since the last change of condition $u$ (in mode $m = 0$) are hidden, while the marker is only observed through noise at each visit date, i.e at each decision date $t_n$ we observe $Y_n = F(X_{t_n}) + \epsilon_n = F(\zeta_{t_n}) + \epsilon_n$ where the $(\epsilon_n)_n$ are independent Gaussian random variables independent from the PDMP and decisions. However, we assume that the death (mode $m = 3$) is perfectly observed.

In this context, the partially observed impulse control problem for our PDMP can be formulated as a POMDP with the following characteristics
<figure><img src="image_26.png" /><figcaption>Figure 24: **Example of observed trajectory in a POMDP framework.** Patient condition modes are now fully hidden, dots indicate (noisy) marker measurements and colors of the dots indicate treatment choices. In this example, the POMDP is controlled by optimization of the value function, and visit dates are allowed to change at each new observation.</figcaption></figure>

* the observation space is $\Omega = \mathbb{R} \times [0, H] \cup \Delta$,

* the observation function does not depend on decision $a$, and has density $o(s' = (x', t'), a, \omega = (y, t)) = f(y - F(x'))1_{t=t'}$ with respect to the Lebesgue measure on $\mathbb{R}$,

* the constraints set is the same as in Example 5.1.3 as constraints did not depend on the hidden part of the process,

* we further assume that the follow-up of the patient starts when they reach their first remission phase, i.e. $b_0 = \delta_{(0,\zeta_0,0)}$.

Figure 24 illustrates an observed trajectory (through noise) of the original PDMP controlled via optimization of the equivalent POMDP value function.

### 5.2.3 Belief states for PDMP

While in the practice of POMDP the belief is often updated through simulations as in the context of *particle filtering* [55, 118], the PDMP context offers a particularly convenient setting for the update of the filters thanks to the Markov property of skeleton chains. Indeed, in the framework of PDMPs, it is common to introduce the filter of the process, and use it for the belief state in the context of POMDPs. General considerations can be found e.g. in [47, 23]. In [30], the authors derive the computations in the specific case of the medical framework example (with slight changes in the model). Here, we briefly expose the steps in the partially observed controlled PDMP framework.

For $0 \le n \le H - 1$, set $\mathcal{F}_n^O = \sigma(Y_k, 0 \le k \le n)$ the $\sigma$-field generated by the observations up to $n$ corresponding to time $t_n$. Let

$$ \Theta_n(B) = \mathbb{P}(X_{t_n} \in B | \mathcal{F}_n^O), $$

for all Borel set $B \subset E_M$, denote the filter or belief process for the hidden part of the process. The standard prediction-correction approach yields a recursive construction for the filter. For any $0 \le n \le H - 1$, conditionally on $Y_{n+1} = y'$, $a = (\ell, r) \in L \times \mathbb{T}$ and $\Theta_n = \theta$, one has $\Theta_{n+1} = \Psi(\theta, y', d)$ with

$$ \Psi(\theta, y', d)(B) = \frac{\int_{E_M} \int_{E_M} f(y' - F(x'))1_B(x') P_r^\ell(dx'|x)\theta(dx)}{\int_{E_M} \int_{E_M} f(y' - F(x'))P_r^\ell(dx'|x)\theta(dx)}, $$
for any Borel subset $B$ of $E_M$. Conditional on being able to compute the integrals, this provides a
theoretical framework for belief updates that does not rely on the quality of simulators. In practice,
computing this filter is not easy and may involve approximations through discretizations.

# 6 Concluding remarks and open questions

We conclude this survey by stating some open questions regarding the control of PDMPs that we believe
can be tackled through the MDP formulation as well as interesting new questions arising in MDP from
this PDMP perspective. As shown in the previous section, the formalism of MDPs makes it possible to
address very challenging problems for controlled PDMPs. We briefly present some of these challenges
in Section6.1. Conversely, deriving MDPs from PDMPs opens new questions for continuous-state space
MDPs and their practical resolution. We briefly present some of these challenges in Section6.2.

## 6.1 Open questions in controlled PDMPs that might be tackled using MDPs

The very well established field of MDPs, together with its wide range of available resolution algorithms,
represents a unique opportunity to tackle open questions that prevail in the framework of impulse control
for PDMPs. A few examples have already been described above, and we provide here a short and
non-exhaustive list of problems that we believe might be formalized and solved in an MDP context.

As mentioned in Sections 2.2 and 5.1.1, the very task of formalizing the definition of an impulse
controlled PDMP in the partially observed framework is rather technical (even in the fully observed
case). Embedding it in the MDP framework provides a very natural environment for this formalization,
since by construction strategies are admissible policies, and value functions are natural elements of MDP
resolutions.

The MDP framework provides a wide range of algorithms for the efficient resolution of controlled PDMPs. As mentioned in Section 5.1 and 5.2, a (partially observed) controlled PDMP can be cast as an MDP with a continuous state-space. From there, several possibilities can be considered: discretizing the initial process to reduce it to a finite state space, for which exact resolution through dynamic programming is possible, or using approximate resolution algorithms designed for MDPs with finite but very large state spaces, or continuous state spaces. Those two strategies have been developed for variants of the medical example, see for instance [30] for an exact resolution of the discretized process in the framework of optimal stopping, [31] for an exact resolution of the discretized process in the general impulse control framework, and [52] for an approximate resolution based on the adaptation of POMCP for PDMPs.

Noting that (partially observed) PDMPs are both easy to simulate and can be cast as POMDPs, one
can take advantage of the huge amount of resources (both in the literature and software) available in the
field of (deep)-Reinforcement Learning (see, e.g. [100, 125]). The field of RL is mature and there exist
several RL application programming interfaces (e.g. Gymnasium [20]), connecting POMDP simulation
models, to libraries of (deep) RL algorithms (e.g. RLlib [97]). The embedding of PDMPs in the MDP
paradigm hence allows to address the scalability issues of standard resolution strategies for PDMPs,
by providing efficient resolution algorithms, including exact and approximate approaches, that provide
excellent performances.

One may also want to consider impulse control for PDMPs where the impulsion does not affect the dy-
namics of the process (for instance deciding how to allocate an observation budget to learn about a process
without affecting trajectories). Here again, while the PDMP framework allows to model continuous-time
trajectories with few parameters that would all be interpretable, the MDP framework should provide the
tools necessary to the formalism and resolution of the control problem of continuous-time trajectories
that are only observed at discrete times.

Finally, the MDP framework also provides the opportunity to solve impulse control problems where
the PDMP model is only partially known, which is different from state partial observability. PDMPs are
typically described by three parameters: the flow $\Phi$, the jump intensity $\lambda$, and the jump kernel $Q$. Solving
MDP problems with ill-known parameters is known as *controlling while learning*, since we have to balance
optimizing the MDP with a currently *wrong* model with *improving our knowledge about the model*. We
have shown in Section 3.3.1 how Bayes-Adaptive MDP allow to solve MDP for which the transition
kernel is unknown or partially known, by modeling the MDP with ill-known parameters as a new MDP
with additional state variables (the hyperparameters). BAMDP are a *control while learning* approach
particularly adapted to solve impulse control problems where the jump kernel $Q$ is only (partially) known.
One can expect that an adaptation of BAMDPs could allow the consideration of other impulse control problems with partially known models, for example, those in which certain parameters of the flow or jump intensity are unknown. Once we have built the BAMDP model, any MDP-solving method (including RL methods can be applied to solve the augmented problem).

## 6.2 Interesting questions in MDPs arising from converted PDMPs

In the MDP/RL community, there is extensive literature on solving problems with finite or continuous time and state spaces (see Section 4). Recent advancements in deep RL have also led to practical solutions for these problems. However, solving MDPs with continuous state spaces or in continuous time remains a significant challenge, and practical algorithms often fail to perform well in these scenarios.

This difficulty is particularly obvious when the transition function, $P(\cdot|s,a)$, is a continuous probability distribution. Monte Carlo algorithms, for instance, generally struggle in such settings because they cannot effectively explore the resulting tree, which has infinitely many successor nodes at each level. Similarly, Reinforcement Learning approaches require approximations of the $Q$-values. While deep RL techniques model $Q$-values using neural networks and are applicable in these cases, they demand immense computational resources and, when relying on real-world data rather than simulations, require vast quantities of data.

Using MDPs derived from piecewise-deterministic Markov processes (PDMPs) offers a way to address complex, realistic problems with models that are both easily parameterized and straightforward to simulate. Additionally, the min-expectation tree resulting from applying Monte Carlo Tree Search (MCTS) to PDMPs maintains a finite, often small branching factor (see e.g. [52]).

In other words, controlled PDMPs serve as an intriguing benchmark for evaluating MDP solution algorithms. Their complexity lies between that of finite MDPs, which are relatively easy to solve, and MDPs with continuous time, state, or decision spaces, where existing solution methods typically lack quality guarantees. Thus, we believe this paper makes a meaningful contribution to the MDP community by addressing these intermediate cases and providing insights into solving such challenging problems.

## References

[1] Anthony Almudevar. A Dynamic Programming Algorithm for the Optimal Control of Piecewise Deterministic Markov Processes. *SIAM Journal on Control and Optimization*, 40(2):525–539, 2001.

[2] Dilip Arumugam and Satinder Singh. Planning to the Information Horizon of BAMDPs via Epistemic State Abstraction. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, *Advances in Neural Information Processing Systems*, volume 35, pages 20482–20497. Curran Associates, Inc., 2022.

[3] C.G. Atkeson and J.C. Santamaria. A comparison of direct and model-based reinforcement learning. In *Proceedings of International Conference on Robotics and Automation*, volume 4, pages 3557–3564 vol.4, April 1997.

[4] Romain Azaïs. A recursive nonparametric estimator for the transition kernel of a piecewise-deterministic Markov process. *ESAIM: Probability and Statistics*, 18:726–749, 2014.

[5] Romain Azaïs and Florian Bouguet, editors. *Statistical Inference for Piecewise-deterministic Markov Processes*. Wiley, 1 edition, August 2018.

[6] Romain Azaïs, François Dufour, and Anne Gégout-Petit. Non-Parametric Estimation of the Conditional Distribution of the Interjumping Times for Piecewise-Deterministic Markov Processes. *Scandinavian Journal of Statistics*, 41(4):950–969, December 2014.

[7] Romain Azaïs and Alexandre Génadot. Semi-parametric inference for the absorption features of a growth-fragmentation model. *TEST*, 24(2):341–360, June 2015.

[8] Romain Azaïs and Alexandre Génadot. A new characterization of the jump rate for piecewise-deterministic Markov processes with discrete transitions. *Communications in Statistics - Theory and Methods*, 47(8):1812–1829, April 2018.
[9] Romain Azaïs and Alexandre Génadot. Estimation of the average number of continuous crossings for non-stationary non-diffusion processes. *Journal of Statistical Planning and Inference*, 198:119-138, January 2019.

[10] Elena Bandini. Optimal control of Piecewise Deterministic Markov Processes: a BSDE representation of the value function. *ESAIM: Control, Optimisation and Calculus of Variations*, 24(1):311–354, 2018.

[11] Vlad Stefan Barbu and Nikolaos Limnios. *Semi-Markov Chains and Hidden Semi-Markov Models toward Applications: Their use in Reliability and DNA Analysis*, volume 191 of Lecture Notes in Statistics. Springer New York, New York, NY, 2008.

[12] J. B. Bardet, A. Christen, A. Guillin, F. Malrieu, and PA. Zitt. Total variation estimates for the TCP process. *Electronic Journal of Probability*, 18:no. 10, 21, 2013.

[13] Camille Baysse, Didier Bihannic, Anne Gégout-Petit, Michel Prenat, Benoîte de Saporta, and Jérôme Sarraco. Maintenance optimization of optronic equipment. *Chemical Engineering Transactions*, 33:709–714, July 2013.

[14] Richard Bellman. Dynamic programming and stochastic control processes. *Inform. And Control (Shenyang)*, 1(3):228–239, September 1958. Publisher: Academic Press.

[15] Michel Benaïm, Stéphane Le Borgne, Florent Malrieu, and Pierre-André Zitt. Quantitative ergodicity for some switched dynamical systems. *Electronic Communications in Probability*, 17(none), January 2012.

[16] D. Bertsekas. Convergence of discretization procedures in dynamic programming. *IEEE Transactions on Automatic Control*, 20(3):415–419, June 1975. Conference Name: IEEE Transactions on Automatic Control.

[17] Joris Bierkens, Alexandre Bouchard-Côté, Arnaud Doucet, Andrew B. Duncan, Paul Fearnhead, Thibaut Lienart, Gareth Roberts, and Sebastian J. Vollmer. Piecewise deterministic Markov processes for scalable Monte Carlo on restricted domains. *Statistics & Probability Letters*, 136:148–154, May 2018.

[18] C. Boutilier, T. Dean, and S. Hanks. Decision-Theoretic Planning: Structural Assumptions and Computational Leverage. *Journal of Artificial Intelligence Research (JAIR)*, 11, July 1999. Publisher: AI Access Foundation.

[19] Adrien Brandejsky, Benoîte de Saporta, and François Dufour. Numerical method for expectations of piecewise deterministic Markov processes. *Communications in Applied Mathematics and Computational Science*, 7(1):63–104, March 2012.

[20] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba. Openai gym. *arXiv preprint arXiv:1606.01540*, 2016.

[21] Cameron B. Browne, Edward Powley, Daniel Whitehouse, Simon M. Lucas, Peter I. Cowling, Philipp Rohlfshagen, Stephen Tavener, Diego Perez, Spyridon Samothrakis, and Simon Colton. A Survey of Monte Carlo Tree Search Methods. *IEEE Transactions on Computational Intelligence and AI in Games*, 4(1):1–43, March 2012. Conference Name: IEEE Transactions on Computational Intelligence and AI in Games.

[22] Nicole Bäuerle and Dirk Lange. Optimal Control of Partially Observable Piecewise Deterministic Markov Processes. *SIAM Journal on Control and Optimization*, 56(2):1441–1462, January 2018.

[23] Nicole Bäuerle and Ulrich Rieder. *Markov Decision Processes with Applications to Finance*. Universitytext. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.

[24] Pablo Samuel Castro and Doina Precup. Using Linear Programming for Bayesian Exploration in Markov Decision Processes. In *International Joint Conference on Artificial Intelligence*, 2007.

[25] Djalil Chafaï, Florent Malrieu, and Katy Paroux. On the long time behavior of the TCP window size process. *Stochastic Processes and their Applications*, 120(8):1518–1534, August 2010.
[26] Caroline Carvalho Chanel, Florent Teichteil-Königsbuch, and Charles Lesire. Multi-target detection and recognition by UAVs using online POMDPs. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 27, pages 1381–1387, 2013.

[27] Hyeong Soo Chang, Jiaqiao Hu, Michael C Fu, and Steven I Marcus. *Simulation-based algorithms for Markov decision processes*, volume 68. Springer, 2013.

[28] Julien Chiquet and Nikolaos Limnios. A method to compute the transition function of a piece-wise deterministic Markov process with application to reliability. *Statistics & Probability Letters*, 78(12):1397–1403, September 2008.

[29] Julien Chiquet, Nikolaos Limnios, and Mohamed Eid. Piecewise deterministic Markov processes applied to fatigue crack growth modelling. *Journal of Statistical Planning and Inference*, 139(5):1657–1667, May 2009.

[30] Alice Cleynen and Benoîte de Saporta. Change-point detection for piecewise deterministic Markov processes. *Automatica*, 97:234–247, November 2018.

[31] Alice Cleynen and Benoîte de Saporta. Numerical method to solve impulse control problems for partially observed piecewise deterministic Markov processes, 2021.

[32] Bertrand Cloez, Benoîte de Saporta, and Maud Joubaud. Optimal stopping for measure-valued piecewise deterministic Markov processes. *Journal of Applied Probability*, 57(2):497–512, June 2020.

[33] Bertrand Cloez, Renaud Dessalles, Alexandre Génadot, Florent Malrieu, Aline Marguet, and Ro main Yvinec. Probabilistic and Piecewise Deterministic models in Biology. *ESAIM: Proceedings and Surveys*, 60:225–245, 2017.

[34] Christiane Cocozza-Thivent. *Processus stochastiques et fiabilité des systèmes*. Number 28 in Mathématiques & applications. Springer-Verl, Berlin Heidelberg Paris [etc.], 1997.

[35] Christiane Cocozza-Thivent. *Markov Renewal and Piecewise Deterministic Processes*, volume 100 of *Probability Theory and Stochastic Modelling*. Springer International Publishing, Cham, 2021.

[36] Manon Costa. A piecewise deterministic model for a prey-predator community. *The Annals of Applied Probability*, 26(6), December 2016.

[37] O. L. V. Costa. Discretizations for the average impulse control of piecewise deterministic processes. *Journal of Applied Probability*, 30(2):405–420, June 1993.

[38] O. L. V. Costa and M. H. A. Davis. Approximations for optimal stopping of a piecewise-deterministic process. *Mathematics of Control, Signals, and Systems*, 1(2):123–146, June 1988.

[39] O. L. V. Costa and M. H. A. Davis. Impulse control of piecewise-deterministic processes. *Mathematics of Control, Signals, and Systems*, 2(3):187–206, September 1989.

[40] O. L. V. Costa and F. Dufour. Stability and Ergodicity of Piecewise Deterministic Markov Processes. *SIAM Journal on Control and Optimization*, 47(2):1053–1077, January 2008.

[41] O. L. V. Costa and F. Dufour. Average Continuous Control of Piecewise Deterministic Markov Processes. *SIAM Journal on Control and Optimization*, 48(7):4262–4291, January 2010.

[42] Oswaldo Luiz Do Valle Costa and François Dufour. *Continuous Average Control of Piecewise Deterministic Markov Processes*. SpringerBriefs in Mathematics. Springer New York, New York, NY, 2013.

[43] Rémi Coulom. Efficient selectivity and backup operators in Monte-Carlo tree search. In *5th International Conference on Computer and Games*, pages 72–83, 2006.

[44] M. H. A. Davis. Piecewise-Deterministic Markov Processes: A General Class of Non-Diffusion Stochastic Models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 46(3):353–376, July 1984.