2.1. Impulse control for hidden PDMPs

Let us first define the class of controlled PDMPs we consider. Let $\mathcal{M}=L\times M$ be a two-dimensional finite set. We will call regimes or modes its elements. We use a product set to distinguish between modes in L that can be controller-chosen and modes in M that cannot. For all regimes $(\ell,m) \in \mathcal{M}$ , let $E_m^\ell$ be a bounded Borel subset of $\mathbb{R}^d$ endowed with a norm $\|\cdot\|$ . Set $E=\{(\ell,m,\mathtt{x}), \ell\in L, m\in M, \mathtt{x} \in E_m^\ell\}$ , and $E^\ell=\{(m,\mathtt{x}), m\in M, \mathtt{x} \in E_m^\ell\}$ for all $\ell\in L$ . A PDMP on the state space E is determined by three local characteristics:

• • The flow $\boldsymbol\Phi(x,t)=(\ell,\boldsymbol\Phi^\ell(m,\mathtt{x},t))=(\ell,m,\Phi_m^\ell(\mathtt{x},t))$ for all $x=(\ell,m,\mathtt{x})$ in E and $t\geq 0$ , where $\Phi_m^\ell\,:\, \mathbb{R}^d \times \mathbb{R}_+ \rightarrow \mathbb{R}^d$ is continuous and satisfies a semi-group property $\Phi_m^\ell(\cdot,t+s)=\Phi_m^\ell(\Phi_m^\ell(\cdot,t),s)$ , for all $t,s\in \mathbb{R}_+$ . This describes the deterministic trajectory between jumps.

• • The jump intensity $\boldsymbol\lambda(x)=\boldsymbol{\lambda}^\ell(m,\mathtt{x})=\lambda_m^\ell(\mathtt{x})$ for all $x=(\ell,m,\mathtt{x})$ in E, where $\lambda_m^\ell\,:\, {E}_m^\ell\rightarrow \mathbb{R}^+$ is a measurable function such that for any $x=(\ell,m,\mathtt{x})$ in E, there exists $\epsilon >0$ such that

  \begin{align*} \int_0^{\epsilon} \lambda_m^\ell(\Phi_m^\ell(\mathtt{x},t)) \, {\textrm{d}} t < +\infty.\end{align*}
  For all $x=(\ell,m,\mathtt{x})$ in E and $t\geq 0$ , set
  \begin{align*}\boldsymbol\Lambda(x,t)&=\boldsymbol\Lambda^\ell(m,\mathtt{x},t)=\Lambda_m^\ell(\mathtt{x},t)=\int_0^{t} \lambda_m^\ell(\Phi_m^\ell(\mathtt{x},s)) \,{\textrm{d}} s.\end{align*}

• • The Markov kernel ${\mathbf{Q}}$ on $(\mathcal{B}({E}),{E})$ represents the transition measure of the process and allows us to select the new location after each jump. It satisfies ${\mathbf{Q}}(\{x\}|x)=0$ for all $x \in {E}$ . We also write $Q_{m}^{\ell}(\!\cdot\,|\mathtt{x})={\mathbf{Q}}(\!\cdot\!|x)={\mathbf{Q}}^\ell(\!\cdot\!|m,\mathtt{x})$ for all $x=(\ell,m,\mathtt{x})\in E$ and add the additional constraint that Q cannot change the value of $\ell$ , as $\ell$ is intended to be controller-chosen, that is, $Q_{m}^{\ell}$ sends $E^\ell$ onto itself.

The formal probabilistic apparatus necessary to precisely define controlled trajectories and to formally state the impulse control problem is rather cumbersome, and will not be used in the sequel. Therefore, for the sake of simplicity, we only present an informal description of the optimal control problem. The interested reader is referred to [Reference Costa and Davis13], [Reference Davis15, section 54], or [Reference Dufour, Horiguchi and Piunovskiy23] for a formal setting. Our optimization problem will be rigorously stated as a POMDP in Section 2.2.

We consider a finite-horizon problem. Let $H>0$ be the optimization horizon, $\delta>0$ a fixed minimal lapse such that $N=H/\delta$ is an integer. Note that $\delta$ is not supposed to be small. A general impulse strategy $\mathcal{S}=(\ell_n,r_n)_{0\leq n\leq N-1}$ is a sequence of non-anticipative E-valued random variables on a measurable space $(\Omega,\mathcal{F})$ and of non-anticipative intervention lapses. In this work, we only consider a subclass of strategies where $\ell_n$ takes values in L and $r_n$ is a multiple of $\delta$ : $r_n\in\mathbb{T}$ , where $\mathbb{T}$ is a subset of $\delta^{1\,:\,N}=\{\delta,2\delta,\ldots, N\delta\}$ that contains $\delta$ . This means that on the one hand the controller can only act on the process by changing its regime, that is, by selecting the local characteristics to be applied until the next intervention, and on the other hand the lapse between consecutive interventions belongs to the finite set $\mathbb{T}$ . The trajectories of the PDMP controlled by strategy $\mathcal{S}$ are constructed recursively between intervention dates $\tau_n$ (defined recursively by $\tau_0=0$ and $\tau_{n+1}=\tau_{n}+r_n$ ) as described in Algorithm 1.In line 5 of Algorithm 1, $S\sim \lambda_m^{\ell}(\mathtt{x})$ means that S has the survival function

\begin{align*}{\mathrm{P}}_x(S>t)={\textrm{e}}^{-\int_0^t \lambda_m^{\ell}(\Phi_m^{\ell}(\mathtt{x},s))\,{\textrm{d}} s}={\textrm{e}}^{-\boldsymbol\Lambda(x,t)}.\end{align*}

Algorithm 1

Simulation of a trajectory of a controlled PDMP between interventions n and $n+1$ from state $x_n=(\ell_n, m_n, \mathtt{x}_n)$

Algorithm 1 Long description

An algorithm for simulating a controlled process trajectory between interventions. The algorithm starts by initializing the state variables s, l, m, and x. It then enters a loop that continues while s is less than r_n. Inside the loop, the algorithm samples a random variable S from a distribution dependent on l and m. It then updates the state variables based on the sampled value of S and other conditions. If s is less than r_n minus s, the algorithm updates l, m, and x based on a function Q. Otherwise, it updates s to be equal to r_n. The loop continues until s is no longer less than r_n.

Note that boundary jumps are not allowed. The strategy $\mathcal{S}$ induces a family of probability measures ${\mathrm{P}}^{\mathcal{S}}_x$ , $x\in E$ , on a suitable probability space $(\Omega,\mathcal{F})$ . Associated to strategy $\mathcal{S}$ , we define the following expected total cost for a process starting at $x\in E$ :

(2.1) \begin{align}\mathcal{J}(\mathcal{S},x)={\mathrm{E}}_x^{\mathcal{S}} \Bigl[ \int_0^{H} c_r(X_s)\,{\textrm{d}} s + \sum_{n=0}^{N-1} c_i({X}_{\tau_n},\ell_n,r_n)+c_t(X_H)\Bigr],\end{align}

where ${\mathrm{E}}_x^{\mathcal{S}}$ is the expectation with respect to ${\mathrm{P}}^{\mathcal{S}}_x$ , $c_r$ is the running cost, $c_i$ the intervention cost, and $c_t$ the terminal cost.

The last ingredient needed to state the optimization problem is to define admissible strategies. Again, the rigorous definition will be given in Section 2.2 in the framework of POMDPs. Informally, decisions can only be taken in view of some discrete-time noisy observations of the process, instead of the exact value of the process at all times. More specifically, we make the following assumptions.

• • Observations are only available at decision times $\tau_n$ .

• • The controller-chosen regimes $\ell\in L$ are observed, while the uncontrolled regimes $m\in M$ are hidden.

• • The Euclidean variable is observed through noise: at time $\tau_n$ , if the state of the process is $x=(\ell,m,\mathtt{x})$ , the controller receives observation $y_n=F(x)+\epsilon_n=F(\mathtt{x})+\epsilon_n$ , where $(\epsilon_n)$ are real-valued independent and identically distributed random variables with density f independent of the controlled PDMP. We further assume that the random variables $y_n$ take values in a compact interval I of the real line.

Denote by $\mathbb{S}$ the set of all admissible strategies. Our aim is to compute an approximation of the value function

\begin{align*}\mathcal{V}(x)=\inf_{\mathcal{S}\in\mathbb{S}}\mathcal{J}(\mathcal{S},x),\end{align*}

and explicitly construct a strategy close to optimality. We also assume that the process starts from a fixed point $x_0=(\ell_0,m_0,\mathtt{x}_0)\in E$ .

2.2. Partially observed Markov decision process

Finding a suitably rigorous way to state an impulse control problem for PDMPs with hidden jumps and under noisy observation is by no means straightforward, especially as regards defining admissible strategies (see, for instance, [Reference Almudevar1] or [Reference Costa and Dufour11, Section 1.1]). This is our first main contribution in this paper. As decision dates are discrete, we use the framework of POMDPs to rigorously state our control problem. In the sequel, we drop the regime $\ell\in L$ from the state x and include it in the action instead to better fit the standard POMDP notation. Denote $E_m=\bigcup_{\ell\in L}E_m^\ell$ and $E_M=\{(m,\mathtt{x}), m\in M, \mathtt{x} \in E_m\}$ . Define on $E_M$ the following hybrid distance, for $x=(m,\mathtt{x})$ and $x'=(m',\mathtt{x}')\in E_M$ :

\begin{align*}\|x-x'\| = (\|{\mathtt{x}}-{\mathtt{x}}'\|) \unicode{x1D7D9}_{m=m'}+\infty \unicode{x1D7D9}_{m \neq m'}.\end{align*}

Let $(\mathbb X,\mathbb A, \mathbb K, R,c,C)$ be the POMDP with the following characteristics.

$\bullet$ The state space is $\mathbb X = \{\xi=(x,y,w)\in (E_M\times \mathbb{O})\}\cup \{\Delta\}$ , with $\mathbb{O}=\{\gamma \in I\times\delta^{1:N}\}$ the observation space. It gathers the values $X_n$ of the hidden PDMP at dates $\tau_n$ as well as the observation process $Y_n$ , and one additional observed variable $W_n$ that is the time elapsed since the beginning. State $\Delta$ is a cemetery state where the process is sent after the horizon H is reached.

$\bullet$ The action space is $\mathbb A = (L\times\mathbb T)\cup \{\check{d}\}$ , where $\check{d}$ is an empty decision that is taken when the horizon H is reached. This purely technical decision sends the process to the cemetery state $\Delta$ .

$\bullet$ The constraints set $\mathbb K\subset \mathbb X\times \mathbb A$ is such that its sections $\mathbb K(\xi)=\{d\in \mathbb A; (\xi,d)\in\mathbb K\}$ satisfy $\mathbb K(\Delta)=\{\check{d}\}$ and $\mathbb K(x,y,w)=\mathbb{K}(w)$ for $(x,y,w)\in\mathbb{X}-\{\Delta\}$ as constraints concern the remaining control time. More specifically, one has:

1. – $\mathbb{K}(w)\subset L\times(\mathbb{T} \cap \{\delta,2\delta,\dots,N\delta-w\})$ if $w\leq \delta(N-1)$ , to force the last intervention to occur exactly at the horizon time H;

2. – $\mathbb{K}(N\delta)=\check{d}$ , denoting that no intervention is possible after the horizon has been reached.

Note that as $\delta\in\mathbb{T}$ , the set $\mathbb{T} \cap \{\delta,2\delta,\dots,N\delta-w\}$ is never empty.

$\bullet$ The controlled transition kernels R are defined as follows: for any bounded measurable function g on $\mathbb{X}$ , any $\xi\in\mathbb{X}$ and $d\in\mathbb{K}(\xi)$ , one has

\begin{align*}Rg(\xi,\check{d})&=g(\Delta),\\Rg(\xi,d)&=\int_{I}\int_{E} g(x',y',w+r)f(y'-F(x'))P({\textrm{d}} x'|x,d)\, {\textrm{d}} y',\end{align*}

if $\xi=(x,y,w)$ , $d=(\ell,r)\neq\check{d}$ , and $Rg(\xi,d)=0$ otherwise, where $P(\!\cdot\!| x,d)$ is the distribution of $X_r$ conditionally on $X_0=x$ under regime $\ell$ , if $d=(\ell,r)$ . Note that since, conditionally on $X_n$ , $Y_n$ is independent of any past $X_k$ and $Y_k$ , kernel R does not depend on y.

$\bullet$ The terminal cost function $C\,:\,\mathbb X\to \mathbb R_+$ satisfies $C(\Delta)=0$ and $C(x,y,w)=c_t(x)$ for $(x,y,w)\in \mathbb{X}-\{\Delta\}$ .

$\bullet$ The cost-per-stage function $c\,:\,\mathbb K\times\mathbb X\to \mathbb R_+$ satisfies $c(\Delta, \check{d},\cdot)=0$ , $c(\xi,\check{d},\cdot)=C(\xi)$ if $\xi\neq \Delta$ . Ideally, to match equation (2.1), for $d=(\ell,r)$ one should choose c as

\begin{align*}c(x,d)={\mathrm{E}}_x^{\ell} \left[ \int_0^{r} c_r(X_s) \,{\textrm{d}} s\right] + c_i(x,d).\end{align*}

However, the integral part has no simple analytical expression and we choose some simpler proxies instead (see Section 5.1).

$\bullet$ The optimization horizon is finite and equals N.

Classically, the sets of observable histories are defined recursively by $H_0=\mathbb{O}\cup\{\Delta\}$ and $H_n=H_{n-1}\times\mathbb{A}\times (\mathbb{O}\cup\{\Delta\})$ . A decision rule at time n is a measurable mapping $g_n\,:\,H_n\rightarrow\mathbb{A}$ such that $g_n(h_n)\in\mathbb{K}(\gamma_n)$ for all histories $h_n=(\gamma_0,d_0,\gamma_1,d_1,\ldots,\gamma_n)$ . A sequence $\pi=(g_n)_{0:N-1}=(g_0,\ldots,g_{n},\ldots, g_{N-1})$ where $g_k$ is a decision rule at time k is called an admissible policy. Let $\Pi_N$ denote the set of all admissible policies. The controlled trajectory of the POMDP following policy $\pi=(g_n)_{0:N-1}\in\Pi_N$ is defined recursively by $\mathbb{X}i_0\in \mathbb{X}$ and, for $0\leq n\leq N-1$ :

• • $A_{n}=g_n(\mathbb{X}i_n)$ ;

• • $\mathbb{X}i_{n+1}\sim R(\cdot | \mathbb{X}i_{n},A_{n})$ .

Note that the cemetery state $\Delta$ ensures that all trajectories have the same length N, even if they do not have the same number of actual decisions ( $d\neq \check{d}$ ). Then one can define the total expected cost of policy $\pi\in\Pi_N$ starting at $\xi_0\in \mathbb{X}$ as

\begin{equation*} J_{\pi}(\xi_0)={\mathrm{E}}^{\pi}_{\xi_0}\left[\sum_{n=0}^{N-1} c(\mathbb{X}i_n,A_n,\mathbb{X}i_{n+1})+C(\mathbb{X}i_N)\right],\end{equation*}

and our control problem corresponds to the optimization problem

\begin{equation*}V(\xi_0)=\inf_{\pi\in \Pi_N}J_\pi(\xi_0).\end{equation*}

Now the problem is rigorously stated, our next aim is to compute an approximation of the value function V and explicitly construct a strategy close to optimality.

3.1. Regularity assumptions

The first set of assumptions will ensure the regularity of kernel P. First, we require the local characteristics of the PDMP to be sufficiently regular. This is a classical assumption, which usually ensures the regularity of the kernel of the inter-jump-time-post-jump-location Markov chain embedded in the PDMP (see, for instance, [Reference de Saporta, Dufour and Zhang18, Section 1.4]). However, here the kernel P is more complex as any number of jumps may occur between two interventions. Therefore, we need an additional restrictive assumption (Assumption 3.5) to limit the number of jumps between consecutive interventions. Denote by $\overline\delta=\max\{r\in \mathbb{T}\}$ the maximum time between consecutive interventions.

Assumption 3.1. The flow $\boldsymbol \Phi$ is Lipschitz continuous: there exists a positive constant $[\Phi]$ such that for all $\ell \in L$ , $x,x'\in E$ , $0\leq t, t'\leq \overline\delta$ , one has

\begin{align*}|\boldsymbol \Phi^\ell(x,t)-\boldsymbol\Phi^\ell(x',t')|&\leq {[\Phi](\|x-x'\|+|t-t'|)}.\end{align*}

Assumption 3.2. The jump intensity $\boldsymbol \lambda$ is Lipschitz continuous along the flow: there exists a positive constant $[\lambda]$ such that for all $\ell\in L$ , x and $x'\in E_M$ , $0\leq t\leq \overline\delta$ , one has

\begin{align*}|\boldsymbol \lambda^\ell(\Phi^\ell(x,t))-\boldsymbol \lambda^\ell(\Phi^\ell(x',t))|\leq [\lambda]\|x-x'\|.\end{align*}

Assumption 3.3. The jump intensity $\boldsymbol \lambda$ is bounded: there exists a positive constant $\|\lambda\|$ such that for all $x\in E_M$ , and $\ell\in L$ , one has

\begin{align*} |\boldsymbol \lambda^\ell(x)|\leq \|\lambda\|.\end{align*}

In the sequel, for all bounded measurable functions h we note that

\begin{align*}{\mathbf{Q}}^\ell h(x)=\int_{E_M}h(x'){\mathbf{Q}}^\ell({\textrm{d}} x',x).\end{align*}

Assumption 3.4. The Markov kernel ${\mathbf{Q}}$ is Lipschitz continuous along the flow in the following sense: there exists a positive constant [Q] such that for all functions $h\in BL_1(E)$ , for all $\ell\in L$ , for all $x,x'\in E_M$ , one has

\begin{equation*}\left|{\mathbf{Q}}^\ell h(x)-{\mathbf{Q}}^\ell h(x')\right|\leq [Q]\|x-x'\|.\end{equation*}

This last assumption is further discussed in Appendix B.2. We also require additional regularity assumptions on the observation process and the cost functions.

Assumption 3.6. There exist non-negative real constants $L_Y$ , $\underline{f}$ , and $\overline{f}$ such that for all $(x,x') \in E_M$ and $y\in I$ , one has

\begin{align*} \big|\,f(y-F({x}))-f(y-F({x'}))\big| &\leq L_Y\|{x-x'}\|,\\0 < \underline{f} \leq f(y-F({x})) &\leq \overline{f}<+\infty.\end{align*}

We also set $L_f=L_Y|I|$ and $B_f=\overline{f}|I|$ , where $|I|=\int_I dy$ is the length of interval I.

Assumption 3.7. There exist non-negative real constants $L_c$ , $L_C$ , $B_c$ , and $B_C$ such that for all x, x^′,x^′′ in $E_M$ and $d\in L\times \mathbb{T}$ , one has

\begin{align*}|C(x)|&\leq B_C,\\|c(x,d,x')|&\leq B_c,\\|C(x)-C(x')|&\leq L_C\|x-x'\|,\\|c(x,d,x')-c(x'',d,x')|&\leq L_c\|x-x''\|,\\|c(x,d,x')-c(x,d,x'')|&\leq L_c\|x'-x''\|.\end{align*}

3.2. Equivalent MDP on the belief state

To convert a POMDP into an equivalent fully observed MDP is classical, therefore details are omitted. The interested reader may consult, for instance, [Reference Bäuerle and Rieder3, Section 5.3] or [Reference Brandejsky, de Saporta and Dufour6, Reference Cleynen and de Saporta8] for similar derivations. For $n\leq N$ , set $\mathcal{F}_n^O=\sigma(Y_k,W_k,0\leq k\leq n)$ , the $\sigma$ -field generated by the observations up to n. Let

\begin{align*}\Theta_n(A)&={\mathrm{P}}(X_{n}\in A | \mathcal{F}_n^O) = {\mathrm{P}}(\mathbb{X}i_{n}\in A\times\{(Y_n,W_n)\} | \mathcal{F}_n^O)\end{align*}

denote the filter or belief process for the unobserved part of the process, for any Borel subset A of $E_M$ . The standard prediction-correction approach yields a recursive construction for the filter.

Proposition 3.1. For any $n\ge 0$ , conditionally on $(Y_{n+1},W_{n+1})=(y',w')$ , $d=(\ell,r)\in L\times\mathbb T$ , and $\Theta_{n}=\theta$ , one has $\Theta_{n+1}=\Psi(\theta,y',w',d)$ with

\begin{align*} \Psi(\theta,y',w',d)(A) &= \frac{\int_{E_M} \int_{E_M}f(y'-F(x'))\unicode{x1D7D9}_A(x')P({\textrm{d}} x'|x,d)\theta({\textrm{d}} x)}{\int_{E_M} \int_{E_M} f(y'-F(x'))P({\textrm{d}} x'|x,d)\theta({\textrm{d}} x)},\end{align*}

for any Borel subset A of $E_M$ .

Let $\mathcal{P}(A)$ denote the set of probability measures on set A. The equivalent fully observed MDP $(\mathbb X',\mathbb A, \mathbb K', R',c',C')$ is defined as follows.

$\bullet$ The state space is a subset $\mathbb{X}'$ of $(\mathcal{P}(E_M)\times \mathbb{O})\cup\{\Delta\}$ satisfying the following constraints: all $\xi=(\theta,y,w)\in\mathbb{X}'$ satisfy $\theta(E_{m_0})\geq ({{\underline{f}}}{{\overline{f}}}^{-1})^{\frac{w}{\delta}\vee 1}{\textrm{e}}^{-w\|\lambda\|}$ , where $m_0$ is the mode of the fixed starting point. It is necessary to restrict the state space to ensure the regularity of our operators (see Appendix B.3). This constraint is technical, and is guaranteed as soon as the process starts in mode $m_0$ (see Appendix A.1).

$\bullet$ The action space is still $\mathbb A = (L\times\mathbb T)\cup \{\check{d}\}$ .

$\bullet$ The constraints set $\mathbb K'\subset \mathbb X'\times \mathbb A$ is such that its sections $\mathbb K'(\xi)=\{d\in \mathbb A;\ (\xi,d)\in\mathbb K'\}$ satisfy $\mathbb K'(\Delta)=\check{d}$ and $\mathbb K'(\theta,y,w)=\mathbb{K}'(w)=\mathbb{K}(w)$ for $(\theta,y,w)\in \mathcal P(E_M)\times\mathbb{O}$ .

$\bullet$ The controlled transition kernels R ^′ are defined as follows: for any bounded measurable function g on $\mathbb{X}'$ , any $\xi\in\mathbb{X}'$ , and $d\in\mathbb{K}'(\xi)$ , one has $R'g(\xi,d)=g(\Delta)$ if $d=\check{d}$ , and

\begin{align*}R'g(\xi,d)&=\int_{E_M}\int_{I}\int_{E_M} g\big(\Psi(\theta,y',w+r,d),y',w+r\big) f(y'\!-F(x'))P({\textrm{d}} x'|x,d)\,{\textrm{d}} y'\theta({\textrm{d}} x) \end{align*}

if $\xi=(\theta,y,w)$ , $d=(\ell,r)\neq \check{d}$ . See Appendix A.1 for the proof that R ^′ maps $\mathbb{X}'$ onto itself.

$\bullet$ The non-negative cost-per-stage function $c'\,:\,\mathbb K'\to \mathbb R_+$ and the terminal cost function $C'\,:\,\mathbb X'\to \mathbb R_+$ are defined by $C'(\Delta)=c'(\Delta, \check{d},\cdot)=0$ , and for $\xi=(\theta,\gamma)\in \mathcal P(E_M)\times \mathbb{O}$ and $d\in \mathbb K'(\xi)$ ,

\begin{eqnarray*} c'(\xi,d)=c'(\theta,\gamma,d)\ &=&\int_{E_M^2} c(x,d,x') P({\textrm{d}} x'|x,d)\theta({\textrm{d}} x), \\ C'(\xi)=C'(\theta,\gamma)&\ =& \int_{E_M} C(x) \theta({\textrm{d}} x).\end{eqnarray*}

$\bullet$ The optimization horizon is still N.

Denote by $(\mathbb{X}i'_n)$ a trajectory of the fully observed MDP. The cost of strategy $\pi\in\Pi_N$ is

\begin{equation*} J'(\pi,\xi'_0) = {{\mathrm{E}}}^{\pi}_{\xi'_0}\left[\sum_{n=0}^{N-1} c'(\mathbb{X}i'_{n},A_n)+C'(\mathbb{X}i'_{N})\right],\end{equation*}

and the value function of the fully observed problem is

\begin{equation*}V'(\xi'_0)=\inf_{\pi\in \Pi_N}J'(\pi,\xi'_0).\end{equation*}

If $\xi'_0=(\delta_{x_0},y_0,w_0)$ , one has $V'(\xi'_0)=V(x_0,y_0,w_0)$ , so that solving the partially observed MDP is equivalent to solving the fully observed one. In addition, the value function V’ satisfies the well-known dynamic programming equations (see, for instance, [Reference Bäuerle and Lange2, Section 5.3]).

Theorem 3.1. For $\xi\in\mathbb{X}'$ , set $v_N'(\xi)= C'(\xi)$ , and for $0\leq n\leq N-1$ , define by backwards induction

\begin{equation*} v_n'(\xi) = \min_{d\in\mathbb{K'}(\xi)} \left\{ c'(\xi,d) + R' v'_{n+1}(\xi,d) \right\}. \end{equation*}

Let $\xi_0'=(\delta_{(m_0,\mathtt{x}_0)},y_0,0) \in \mathbb{X}'$ . Then we have $v'_0(\xi_0') = V'(\xi_0')=V((m_0,\mathtt{x}_0),y_0,0).$

3.3. Approximations of the value function

Our approximation of the value function is based on discretizations of the underlying state spaces.

3.3.1. First discretization.

Let $\Omega=\{\omega^1,\ldots,\omega^{n_{\Omega}}\}$ be a finite grid on $E_M$ containing at least one point in each mode $m\in M$ . Let $p_\Omega$ denote the nearest-neighbour projection from $E_M$ onto $\Omega$ for the distance defined in Section 2.1. In particular, $p_\Omega$ preserves the mode. Set $\Omega_m=E_m\cap\Omega$ for all $n\in M$ . Let $(C_i)_{1\leq i\leq K}$ be a Voronoi tessellation of $E_M$ associated to $\Omega$ . Namely, $(C_i)_{1\leq i\leq n_{\Omega}}$ is a partition of $E_M$ such that for all $1\leq i\leq n_{\Omega}$ , one has

\begin{equation*} C_i\subset \{x\in E_M\, ;\, \|x-\omega^i\|\leq \|x-\omega^j\| \ \forall j\in 1\,:\,n_\Omega\}.\end{equation*}

We will denote by $\mathcal{D}_i$ the diameter of cell $C_i$ : $\mathcal{D}_i=\sup \{\|x-x'\| ;\, x,x' \in C_i \}$ . As stated in Section 2.1, the state spaces are supposed to be bounded. Denote $\|E_M\|=\sup\{\|x\|, x\in E_M\}$ . Hence, all $\mathcal{D}_i$ are finite and bounded by $\|E_M\|$ .

We define the controlled kernels $\bar P$ from $E_M\times (L\times \mathbb{T})$ onto $\Omega$ as

\begin{equation*} \bar P(\omega^j | x,d) =P(C_j|x,d),\end{equation*}

for all $x\in E$ , $d\in L\times \mathbb{T}$ , and $1\leq j\leq n_{\Omega}$ . In particular, the restriction of $\bar P$ to $\Omega$ is a controlled Markov kernel on $\Omega$ . We now replace kernel P by kernel $\bar P$ in the dynamic programming equations of Theorem 3.1 in order to obtain our first approximation.

Let $\bar{\mathbb{X}}'\subset\{\xi=(\theta,\gamma)\in \mathbb{X}'\,;\, \theta(\Omega)=1, \}\cup\{\Delta\}$ satisfying the following constraints: all $\bar\xi=(\bar\theta,y,w)\in\bar{\mathbb{X}}'$ satisfy $\bar\theta(\Omega_{m_0})\geq ({{\underline{f}}}{{\overline{f}}}^{-1})^{\frac{w}{\delta}\vee 1}{\textrm{e}}^{-w\|\lambda\|}$ , where $m_0$ is the mode of the fixed starting point. Let $\bar\Psi$ be the approximate filter operator defined by replacing the integrals with respect to P in the filter operator $\Psi$ defined in Proposition 3.1 by integrals with respect to $\bar P$ . The approximated controlled transition kernels $\bar R'$ are defined as follows: for any bounded measurable function g on $\bar{\mathbb{X}}'$ , any $\xi\in\bar{\mathbb{X}}'$ and $d\in\mathbb{K}'(\xi)$ , one has $\bar R'g(\xi,d)=g(\Delta)$ if $d=\check{d}$ , and

\begin{align*}\bar R'g(\xi,d)&=\int_{I}\sum_{\omega^i\in\Omega}\sum_{\omega^j\in\Omega} g\big(\bar{\Psi}(\bar{\theta},y',w+r,d),y',w+r\big)f(y'-F(\omega^j))\bar P(\omega^j|\omega^i,d)\bar{\theta}(\omega^i)\,{\textrm{d}} y'\end{align*}

if $\xi=(\bar{\theta},y,w)$ , $d=(\ell,r)$ . Note that this kernel does not depend on y, and see Appendix A.2 for the proof that $\bar R$ sends $\bar{\mathbb{X}}'$ onto itself. The approximated cost-per-stage function is defined by $\bar c'(\Delta,\check{d})=0$ and

\begin{align*} \bar c'(\bar{\theta},\gamma,d)&=\sum_{\omega^i\in\Omega} \sum_{\omega^j\in\Omega}c(\omega^i,d,\omega^j) \bar P(\omega^j|\omega^i,d)\bar{\theta}(\omega^i),\end{align*}

for all $\xi=(\bar{\theta},\gamma)\in\bar{\mathbb{X}}'$ and $d\in\mathbb{K}'(\xi)$ . Finally, for all $\xi\in\bar{\mathbb{X}}'$ , set $\bar v_N'(\xi)= \ C'(\xi)$ , and for $0\leq n\leq N-1$ , define by backwards induction

\begin{eqnarray*} \bar v_n'(\xi)&=&\min_{d\in\mathbb{K}(\xi)} \left\{ \bar c'(\xi,d) + \bar R' \bar v'_{n+1}(\xi,d) \right\}. \end{eqnarray*}

If the grid $\Omega$ is precise enough, $\bar P$ should be a good approximation of kernel P and thus one can expect that functions $\bar v'_n$ are good approximations of our value functions $v'_n$ . Indeed, we have the following result.

Theorem 3.2. Under Assumptions 3.1–3.7, for all $\xi'\in\bar{\mathbb{X}}'$ , one has

\begin{align*} |v'_N(\xi') - \bar{v}_N'(\xi')| & = 0, \\ |v'_n(\xi')- \bar{v}_n'(\xi')| &\leq C_{v_n'}\sup_{j \in \{1,\dots,n_\Omega\}} \mathcal{D}_j, \quad 0\leq n < N,\end{align*}

where $C_{v'_n}$ depends only on n, N, $\delta$ , and the regularity constants of the parameters.

The proof of this result is given in Appendix B.7 and is based on the form of the kernels, their regularity (especially that of kernel P), and the dynamic programming equation.

The main gain with this first approximation is that the filter operator $\bar \Psi$ now involves only finite weighted sums and therefore the corresponding approximate filter process can now be simulated. However, functions $\bar v'_n$ still cannot be computed because of the continuous integration in y’ and because the state space $\bar{\mathbb{X}}'$ is still continuous. We now proceed to a second discretization in order to obtain numerically tractable approximations of our value functions.

3.3.2. Second discretization.

Let $\Gamma=\{\rho^1,\ldots,\rho^{n_{\Gamma}}\}$ be a finite grid on $\bar{\mathbb{X}}'$ containing at least $\Delta$ , one element $\xi=(\theta,\gamma)$ such that $\theta(E_M)=1$ . Let $p_\Gamma$ denote the nearest-neighbour projection from $\bar{\mathbb{X}}'$ onto $\Gamma$ . Here, the distance on $\bar{\mathbb{X}}'$ is such that $d(\xi, \xi')=d_\Omega(\bar{\theta},\bar{\theta}')+|y-y'|$ if $\xi=(\bar{\theta},y,w),\ \xi'=(\bar{\theta}',y',w')$ , with $w=w'$ , $d(\xi,\xi')= +\infty$ otherwise, where $d_\Omega$ is defined as

\begin{align*}d_\Omega(\bar{\theta},\bar{\theta}')&=\left[\sum_{k=1}^{n_{\Omega}}|\bar{\theta}(\omega^k)-\bar{\theta}'(\omega^k)|^2\right]^{1/2},\end{align*}

meaning that $\mathcal P(\Omega)$ is identified with the ${n_{\Omega}}$ -dimensional simplex endowed with the $L^2$ -norm. Let $(\bar C_i)_{1\leq i\leq {n_{\Gamma}}}$ be a Voronoi tessellation of $\bar{\mathbb{X}}'$ associated to the grid $\Gamma$ . We will denote by $\mathcal{\bar D}_i$ the diameter of cell $\bar C_i$ . Note that the space $\bar{\mathbb{X}}'$ is compact, so that all $\mathcal{\bar D}_i$ are finite and bounded by $1+|I|$ . The controlled kernels $\hat R'$ are defined on $\bar{\mathbb{X}}'$ as follows: for any $\rho^{\,j}\in\Gamma$ , any $\xi\in\bar{\mathbb{X}}'$ and $d\in\mathbb{K}'(\xi)$ , one has

\begin{equation*} \hat R'(\rho^{\,j} | \xi,d) = \bar R'(\bar C_j|\xi,d).\end{equation*}

In particular, the restriction of $\hat R'$ to $\Gamma$ is a family of controlled Markov kernels on $\Gamma$ . For $1\leq j\leq {n_{\Gamma}}$ , set $ \hat v_N'(\rho^{\,j})= C'(\rho^{\,j})$ , and for $0\leq n\leq N-1$ and $1\leq j\leq {n_{\Gamma}}$ , define by backwards induction

\begin{eqnarray*} \hat v_n'(\rho^{\,j})&=&\min_{d\in\mathbb{K}'(\rho^{\,j})} \left\{ \bar c'(\rho^{\,j},d) + \hat R' \hat v'_{n+1}(\rho^{\,j},d) \right\}. \end{eqnarray*}

Again, if the grid $\Gamma$ is precise enough, $\hat R'$ should be a good approximation of kernel $\bar R'$ and thus one can expect that functions $\hat v'_n$ are good approximations of functions $\bar v'_n$ and thus of our value functions $v'_n$ . More precisely, we have the following result.

Theorem 3.3. Under Assumptions 3.6 and 3.7, for all points $\rho^{\,j}\in \Gamma$ , one has

\begin{align*} |\hat v'_N(\rho^{\,j}) - \bar{v}_N'(\rho^{\,j})| & =0, \\ |\hat v'_n(\rho^{\,j})- \bar{v}_n'(\rho^{\,j})| &\leq C_{\bar v_n'}\sup_{j \in \{1,\dots,{n_{\Gamma}}\}} \mathcal{\bar D}_j,\end{align*}

where $C_{\bar v'_n}$ depends only on n, N, $\delta$ , and the regularity constants of the parameters.

The proof of this result is given in Appendix C.3. The main gain with this second approximation is that integration against kernel $\hat R'$ boils down to computing a finite weighted sum, and functions $\hat v'_n$ are defined on a finite state space which makes the dynamic programming equations fully tractable numerically.

3.4. Candidate strategy

We can now construct a computable strategy using the fully discretized value function. The idea is to first build an approximate filter using the operator $\bar \Psi$ and project the resulting filter together with the current observation onto grid $\Gamma$ . Then one selects the next decision using the dynamic programming equation on the grid $\Gamma$ . More precisely, suppose that the process starts from point $\xi_0=(m_0,\mathtt{x}_0,y_0,0)$ , such that the initial observation is $\gamma_0=(y_0,0)$ . One can recursively compute an approximate filter $(\bar\theta_n)$ and the corresponding decisions $(d_n)$ as follows. First, set

\begin{align*}\bar\theta_0=\delta_{(m_0,\mathtt{x}_0)},\quad d_0(\bar \theta_0,\gamma_0)=\mathrm{argmin}_{d\in\mathbb{K}'(\bar\theta_0,\gamma_0)} \big\{\bar c'(p_{\Gamma}(\bar \theta_0,\gamma_0),d)+\hat R'\hat v'_{1}(p_{\Gamma}(\bar \theta_0,\gamma_0),d) \big\}.\end{align*}

Suppose one has constructed the sequence $(\bar\theta_n, d_n)$ up to stage $k-1$ . Then after receiving the kth observation $\gamma_k$ , the next approximate filter and decisions are

\begin{align*}\bar\theta_k=\bar\Psi(\bar\theta_{k-1},\gamma_k,d_{k-1}),\ d_k(\bar \theta_k,\gamma_k)=\mathrm{argmin}_{d\in\mathbb{K}'(\bar\theta_0,\gamma_0)} \big\{\bar c'(p_{\Gamma}(\bar \theta_k,\gamma_k),d)+\hat R'\hat v'_{1}(p_{\Gamma}(\bar \theta_k,\gamma_k),d) \big\}.\end{align*}

with the convention that the last (not required) decision is $d_N=\check{d}$ . Note that it should be better to use operator $\bar\Psi$ on $\bar\theta_{k-1}$ then project the result onto grid $\Gamma$ than to use operator $\hat \Psi$ (obtained by replacing integration with respect to $\theta$ by integration with respect to $p_\Gamma(\theta)$ in the definition of $\bar\Psi$ ) as it should generate a smaller error. Although it is reasonable to think that this strategy should be close to optimality, it is an open problem to actually prove it as the sequence $(\bar\theta_n,\gamma_n)$ is not generated by the kernel $\bar R'$ . Indeed, we use here the observations generated by the original sequence $X_n$ with kernel P and not that generated by kernel $\bar P$ , as in practice only the original observations are available. Its performance is assessed in Section 5 on a specific example.

4.1. Optimization of patient follow-up

We consider an example of patient follow-up that we will use in the simulation study in Section 5. The mode m corresponds to the overall state of the patient ( $m=0$ , sound; $m=1$ , disease 1; $m=2$ , disease 2; $m=3$ : death of the patient). The Euclidean variable is two-dimensional $\mathtt{x}=(\zeta,u)$ . Variable $\zeta\in[K_i,K_s]$ corresponds to some marker of the disease that can be measured, and $u\geq 0$ the time since the last jump in order to encompass semi-Markov dynamics. The control $\ell$ corresponds to the medical treatment ( $\ell=\emptyset$ , no treatment; $\ell=a$ , efficient for disease 1 and slows the progression of disease 2; $\ell=b$ , efficient for disease 2 and slows the progression of disease 1). Decision dates correspond to visits to the medical centre when the marker is measured and a new treatment is selected and applied until the next visit. Horizon H is 2400 days with possible visits every 15, 30, or 60 days ( $\delta=15$ , $\overline\delta=60$ ).

If treatment $\emptyset$ is applied, the patient may randomly jump from $m=0$ to any of the two disease states $m\in\{1,2\}$ . In the sound state $m=0$ , the marker level is constant. In the disease states ( $m=1$ or $m=2$ ), the marker level follows an increasing flow and no other jump is possible.

If treatment a is applied, the patient may only randomly jump from $m=0$ to the other disease state $m=2$ . In the disease state $m=1$ , the marker follows a decreasing flow and may randomly jump to the other disease state $m=2$ . In the disease state $m=2$ , the marker follows an increasing flow and no other change of state is possible. The effects of treatment b are similar: decrease of the marker in disease $m=2$ , increase in disease $m=1$ . The specific values of the local characteristics of the PDMP and proof of their verifying the assumptions of Section 3.1 can be found in the supplementary material, and the state graph for mode transitions is given in Figure 2.

Figure 2.

State graph for the changes of mode in the medical example. Letters indicate under which treatments the jumps are possible. Black arrows indicate changes to remission phase, red arrows changes to standard relapse, and blue arrows changes to therapeutic escape

4.2. Examples in biology

One may consider other kinds of therapeutic treatments, such as in [Reference Pasin, Dufour, Villain, Zhang and Thiébaut30] where the treatments correspond to interleukin injections in HIV-infected patients under antiretroviral treatment that must be suitably dosed and timed in order to maintain lymphocytes levels over a prescribed limit. In this example, the PDMP has no free modes and a controller-chosen mode corresponding to non-proliferating and proliferating lymphocytes, which dynamics can be influenced by punctual interleukin injections. Natural jumps correspond to the efficiency duration of an injection. Hence, there is at most one jump between two injections. It is natural to suppose that the lymphocytes are observed through noisy proxies at discrete dates.

Other examples in biology may cover, for instance, gene expression or chemostats. The simplest PDMP models for gene expression include a (possibly controller-triggered) promoter that can be on or off, and continuous levels of RNA and proteins (see, for instance, [Reference Herbach, Bonnaffoux, Espinasse and Gandrillon26]). If only the controller can turn the promoter on or off, there are no natural jumps. One may also consider a network of interacting genes with partial control over promoters and feedback loops. Chemostats are bioreactors used to cultivate some micro-organisms. They can be modelled by PDMPs (see, for instance, [Reference Fritsch, Harmand and Campillo24]), and are controlled by tuning the flow of fresh medium. Again, it is natural to consider hidden information in these frameworks.

4.3. Dynamic reliability

The field of dynamic reliability is particularly suitable for controlled PDMP models, as explained in the introduction. The number of jumps in a given time interval is usually bounded for non-reparable systems, for instance as soon as there is a finite number of components, each with a single nominal state and a finite number of successive possible degradation states up to total failure. It is also common not to observe the state of each component at all times.

The simplest examples in this field are pure jump Markov processes. These fit our assumptions as the flow and intensity are constant in every mode, and the state space is discrete, hence all kernels are regular. See [Reference Makis and Jiang28] for such an example in the simpler case of optimal stopping. The next important subclass is when the mode corresponds to the states of the components, and the Euclidean variables to their age or use-time (see, for instance, [Reference Baysse, Bihannic, Gégout-Petit, Prenat, de Saporta and Saracco4, Reference de Saporta, Dufour and Zhang19] or [Reference de Saporta, Dufour and Zhang18, Section 1.8.4] for a repair-workshop model with similar variables). In this case, the control $\ell$ may correspond to a change of settings.

A typical example of a more complex controlled PDMP in this field is the heated hold-up tank, a well-known test case for dynamic reliability (see, for instance, [Reference Devooght21, Reference Marseguerra and Zio29, Reference Zhang, Dufour, Dutuit and Gonzalez35]. The system consists of a tank containing a fluid whose level is controlled by three components: two inlet pumps and one outlet valve. A thermal power source heats up the fluid. The failure rate of the components depends on the temperature, the position of the three components monitors the liquid level in the tank, and in turn, the liquid level determines the temperature. The pumps and valve can be either functioning or stuck, which corresponds to the free modes m, and turned on or off (if functioning), corresponding to the controller-chosen mode $\ell$ . In any setting $\ell$ , there are at most three jumps (corresponding to the failure of the pumps and valve). The flow and jump parameters have analytical expressions and satisfy our regularity assumptions.

5.1. Choice of the cost functions

The candidate policy depends on the cost functions, therefore the latter must be carefully chosen. We denote by $c_i(x,d)=C_V$ a fixed cost per visit that takes into account an emotional burden for the patient and health-care expenses for the check-up. For the counterpart of the integral of the running cost, we choose en expression of the form

\begin{align*}\tilde c(x,d,x')=\kappa |\zeta'-K_i|r+\beta r \unicode{x1D7D9}_{\{\ell\neq\emptyset\}}.\end{align*}

Parameter $\beta>0$ represents the cost of the treatment, $\kappa>0$ is a scale parameter and $ |\zeta'-K_i|r$ is a (crude) proxy of the integral of the process. Our cost function is then $c(x,d,x')=C_V+\tilde c(x,d,x')$ . Calibrating those parameters is a difficult task that is not discussed here.

5.2. Construction of the grids

Because decisions influence the dynamics of the process, grids cannot be constructed with techniques based on simulations, such as quantization. For instance, on a horizon of 2400 days with possible visits every 15, 30, or 60 days, this leads to approximately $10^{152}$ possible strategies. Grids therefore have to be chosen by expert knowledge, and transition probabilities computed accordingly. In the supplementary materials we discuss the construction of parsimonious grids that provide satisfying performance.

5.3. Strategies in competition and performance criteria

To evaluate the performance of our approach (OS in the comparison tables) we compare our work with several other strategies.

The (unachievable) gold standard strategy is the see all (SA) where decisions are taken while observing the process $X_{t_n}$ by choosing the optimal treatment for the current mode. In this SA strategy, we do not allow choice of the next visit date.

The filter strategy corresponds to choosing the optimal treatment for the estimated mode: at each new observation $y_{n}$ the approximate filter $\bar\theta_{n}$ is computed based on the first discretization. Then the most probable mode $\hat m_{n}$ is obtained as

\begin{align*}\hat m_{n}=\arg\max_{m=0:3}\sum_{j; \omega^j=(m,\cdot)}\bar\theta_{n}^j.\end{align*}

Note that to set up this strategy, only the first discretization is needed. As this discretization is the least computationally extensive, it might be worth designing larger grids $\Omega$ to obtain better mode estimates. However, this strategy does not take into account knowledge of the dynamics of the process, hence we expect it to have higher costs. It is not able to select the next decision date either.

The standard strategy is classically used in the hospitals and is based on a threshold $s_{\mathrm{rel}}$ for relapse, and on a fixed duration for treatment tlength. While the observations remain below $s_{\mathrm{rel}}$ , the patient does not receive treatment and visits are scheduled every 2 months. When $s_{\mathrm{rel}}$ is reached, the practitioner gives treatment b (corresponding to the most frequent relapse type 2) and the next visit is scheduled in 15 days. If at the next visit the observed marker level has decreased, treatment b is maintained with visits every 15 days for a total duration of tlength days. Otherwise the practitioner changes treatment for treatment a (with visits every 15 days, and same total duration). Once treatment is stopped, visits are scheduled every 2 months and the strategy is repeated.

The fixed dates optimal strategies are the strategies based on our discretization and dynamic programming approach where the choice of the next visit date is not allowed. We will investigate the fixed dates every 15 and 60 days (FD-15 and FD-60 in the comparison tables).

The performance criterion is the real cost of each strategy evaluated on the real process X (averaged over 1000 simulations).

5.4. Results

First, note that all code is available online at https://github.com/acleynen/PDMP-control. All results presented here are based on 1000 simulations. The exact specifications can be found in the supplementary materials.

A first interesting result relates to the allocation of the computational burden. We compared, for a fixed number of elements of the second grid ( $n_\Gamma$ ), the performance of our approach with different numbers of points in the first grid ( $n_\Omega$ ), the smallest ones being included in the largest. The performance of the approach drastically improves with a higher number of points in $\Omega$ , leading to a very crude discretization of $\mathcal{P}(\Omega)$ , compared to a less refined discretization of $\mathbb{X}$ but leading to a more refined discretization of $\mathcal{P}(\Omega)$ . In some cases, adding elements to $\Gamma$ in fact worsened the results, as novel elements allowed consideration of filters with less certainty in the values of the process, most often leading to worse decisions. This is illustrated in Table 5 of the supplementary materials.

We therefore compared our approach with competing strategies on the largest $\Omega$ grid, containing 133 elements in $\Omega$ and 144 elements in $\Gamma$ . Figure 3 illustrates a controlled trajectory with treatments indicated by shape and disease indicated by colour, for our approach. One can note that in most cases, the first time the projected filter has a most likely element in a disease mode, treatment is not administered but the time lapse to the next visit is reduced. Only on the next visit will treatment be allocated. This explains the differences of performance between the different strategies given in Table 1 and shows the added value of taking into account the knowledge of the dynamics of the process. It performs significantly better.

Figure 3.

Example of a controlled (OS) simulated strategy. The top panel illustrates the (hidden) values of the true process, with modes indicated by colour, and allocated treatment indicated by point shape. The middle panel illustrates the observed trajectory (noisy version of the Euclidean variable from the panel above). The bottom panel illustrates the most likely element of $\Omega$ estimated by the filter of the process

Table 1.

Strategy comparisons. Performance of each strategy measured in terms of Monte Carlo costs.

Table 1 Long description

A table comparing the performance of different strategies in terms of visits and real cost. The table has 7 rows and 3 columns. The columns are labeled Visits, Real cost (sd), and Choice. The rows are labeled OS, Filter, See all, and Standard. The values in the table are as follows: Row 1: Choice, 203 (3.03). Row 2: OS 15 days, 291.36 (2.66). Row 3: OS 60 days, 215.98 (3.2). Row 4: OS 15 days, 321.28 (2.2). Row 5: Filter 60 days, 223.8 (4.03). Row 6: Filter 15 days, 203.15 (0.89). Row 7: See all 60 days, 113.72 (1.79). Row 8: Standard, 301.09 (5.69).