Application of the interacting particle system method to piecewise deterministic Markov processes used in reliability

TL;DR

This paper introduces IPS+M, an adaptation of the interacting particle system method, which improves variance reduction in reliability assessments of concentrated PDMPs by preventing repetitive trajectories, outperforming the standard IPS.

Contribution

The paper proposes IPS+M, a modified IPS method that enhances variance reduction for concentrated PDMPs by conditioning propagation to avoid identical trajectories, with proven lower variance.

Findings

IPS+M always yields lower variance estimators than IPS.

Simulation confirms IPS+M's improved efficiency on a two-component system.

Concentrated PDMPs pose challenges for standard IPS, which IPS+M addresses.

Abstract

Variance reduction methods are often needed for the reliability assessment of complex industrial systems, we focus on one variance reduction method in a given context, that is the interacting particle system method (IPS) used on piecewise deterministic Markov processes (PDMP) for reliability assessment . The PDMPs are a very large class of processes which benefit from high modeling capacities, they can model almost any Markovian phenomenon that does not include diffusion. In reliability assessment, the PDMPs modeling industrial systems generally involve low jump rates and jump kernels favoring one safe arrival, we call such model a "concentrated PDMP". Used on such concentrated PDMPs, the IPS is inefficient and does not always provide a variance reduction. Indeed, the efficiency of the IPS method relies on simulating many different trajectories during its propagation steps, but unfortunately concentrated PDMPs are likely to generate the same deterministic trajectories over and over. We propose an adaptation of the IPS method called IPS+M that reduces this phenomenon. The IPS+M consists in modifying the propagation steps of the IPS, by conditioning the propagation to avoid generating the same trajectories multiple times. We prove that, compared to the IPS, the IPS+M method always provides an estimator with a lower variance. We also carry out a quick simulation study on a two-components system that confirms this result.

Tables

Table 1. Table 1: Empirical means and empirical variances on 100 100 100 runs with N = 10 5 𝑁 superscript 10 5 N=10^{5} for the MC, the IPS and the IPS+M methods

			MC	IPS	IPS+M
		$\hat{p}$	$\mathrm{\hspace{0.17em}2.71}\times {10}^{-5}$
		${\hat{\sigma }}^2$	$2.90\times {10}^{-10}$
	$n=5$	$\hat{p}$		$\mathrm{\hspace{0.17em}2.86}\times {10}^{-5}$	$2.70\times {10}^{-5}$
		${\hat{\sigma }}^2$		$\mathrm{\hspace{0.17em}1.78}\times {10}^{-9}$	$1.37\times {10}^{-10}$
	$n=10$	$\hat{p}$		$\mathrm{\hspace{0.17em}2.85}\times {10}^{-5}$	$2.64\times {10}^{-5}$
		${\hat{\sigma }}^2$		$\mathrm{\hspace{0.17em}1.08}\times {10}^{-9}$	$1.07\times {10}^{-10}$

Table 2. Table 2: Empirical means and empirical variances on 50 50 50 runs with N = 10 5 𝑁 superscript 10 5 N=10^{5} for the MC, the IPS and the IPS+M methods

			MC	IPS	IPS+M
		$\hat{p}$	$\mathrm{\hspace{0.17em}1.12}\times {10}^{-4}$
		${\hat{\sigma }}^2$	$1.12\times {10}^{-9}$
	$n=5$	$\hat{p}$		$\mathrm{\hspace{0.17em}1.75}\times {10}^{-4}$	$1.12\times {10}^{-4}$
		${\hat{\sigma }}^2$		$\mathrm{\hspace{0.17em}3.08}\times {10}^{-8}$	$4.37\times {10}^{-9}$

Equations

$$\frac{d{X}_t}{dt}=F_{M_t}(X_t).$$

$$\forall t\in[0,T),\quad Z_{s+t}=\big{(}X_{s+t},M_{s}\big{)}=\big{(}\phi_{(X_{s},M_{s})}(t),M_{s}\big{)}.$$

$$\forall t\in [0,T),Z_{t+s}={\Phi }_{Z_s}(t).$$

$$E=\bigcup _{m\in \mathbb{M}}{E}_m.$$

$$\lambda (Z_{s+t})=\sum _{m^{+}\in \mathbb{M}}{\lambda }_{M_s\to m^{+}}(X_{s+t})$$

$$\mathbb{P}(T\leq t|Z_{s}=z)=\left\{\begin{array}[]{cr}1-\exp\left[-\Lambda_{z}(t)\right]&\mbox{ if }t<t^{*}_{z}\,,\\ 1&\mbox{ if }t\geq t^{*}_{z}\,.\\ \end{array}\right.$$

$$\forall B\in \mathcal{B}({\mathbb{R}}^{+}),{\mu }_z(B)$$

$$\forall B\in \mathcal{B}(E),\mathbb{P}\left(Z_T^{+}\in B|Z_T^{\text{-}}=z^{\text{-}}\right)$$

$$\forall B\in \mathcal{B}(E),{\nu }_{z^{\text{-}}}(B)=\sum _{m^{+}\in \mathbb{M}}{\delta }_{z_a(z^{\text{-}},m^{+})}(B).$$

$$\mathbb{P}\left(Z_T^{+}\in B|Z_T^{\text{-}}=z^{\text{-}}\right)$$

$$\forall z^{\text{-}}\in E,B\in \mathcal{B}(E),{\mathcal{K}}_{z^{\text{-}}}(B)$$

$$\text{mbox}and{K}_{z^{\text{-}}}(z^{+})$$

$$K_{z^{\text{-}}}(z^{+})=\text{mathbbm}{1}_{x^{+}=x_a(z^{\text{-}},m^{+})}\prod _{i\in ask(z^{\text{-}},z^{+})}{\left({\kappa }_i(z^{-})\right)}^{fod(i,m^{\text{-}},m^{+})}{\left(1-{\kappa }_i(z^{-})\right)}^{(1-fod(i,m^{\text{-}},m^{+}))}.$$

$$\displaystyle\mathbb{P}_{z_{o}}\Big{(}\mathbf{Z_{t}}\in\Theta_{t}^{-1}(B)\Big{)}=$$

$$\times d{\delta }_{t_n^{\ast }}(t_n)d{\nu }_{z_n^{-}}(z_n)d{\mu }_{z_{n-1}}(t_{n-1})...d{\nu }_{z_1^{-}}(z_1)d{\mu }_{z_o}(t_0),$$

$$\displaystyle\mathbb{P}\big{(}\mathbf{Z}_{t}=\mathbf{z}_{t}\big{)}=\prod_{k=0}^{n}\exp\Big{[}-\Lambda_{z_{s_{k}}}(t^{*}_{z_{s_{k}}})\Big{]}\prod_{k=1}^{n}K_{z_{s_{k}}^{-}}(z_{s_{k}})>0,$$

$$\forall B\in {\mathcal{S}}_t,{\zeta }_{z_0,t}({\Theta }^{-1}(B))=$$

$$\mathbb{P}_{z_{s_{k}}}\left(T_{k}=t^{*}_{z_{s_{k}}}\right)=\exp\Big{[}-\Lambda_{z_{s_{k}}}(t^{*}_{z_{s_{k}}})\Big{]}\simeq 1.$$

$$K_{z_{s_k}^{-}}(z_{s_k})\simeq 1.$$

$$\displaystyle\mathbb{P}\big{(}\mathbf{Z}_{t}=\mathbf{z}_{t}\big{)}=\prod_{k=0}^{n}\exp\Big{[}-\Lambda_{z_{s_{k}}}(t^{*}_{z_{s_{k}}})\Big{]}\prod_{k=1}^{n}K_{z_{s_{k}}^{-}}(z_{s_{k}})\simeq 1.$$

$$p_h=\mathbb{E}[h({\mathbf{Z}}_{t_f})].$$

$$p_{\mathscr{D}}=\mathbb{P}_{z_{0}}\big{(}\mathbf{Z}_{t_{f}}\in\mathscr{D}\big{)}=\mathbb{E}_{z_{0}}\big{[}\mathbbm{1}_{\mathscr{D}}(\mathbf{Z}_{t_{f}})\big{]}.$$

$$\forall {\mathbf{z}}_{{\tau }_k}\in {\mathbf{E}}_{{\tau }_k},B\in {\mathcal{E}}_{{\tau }_{k+1}},\mathbb{P}\left({\mathbf{Z}}_{{\tau }_{k+1}}\in B|{\mathbf{Z}}_{{\tau }_k}={\mathbf{z}}_{{\tau }_k}\right)=V_k(B|{\mathbf{z}}_{{\tau }_k}).$$

$$\forall {\mathbf{z}}_k\in {\mathbf{E}}_k,G_k({\mathbf{z}}_k)$$

$${\tilde{\eta }}_k(d{\mathbf{z}}_{{\tau }_k})$$

$$\forall B\in {\mathcal{E}}_{{\mathbf{\tau }}_{\mathbf{k}}},{\tilde{\eta }}_k(B)$$

$${\eta }_{k+1}(d{\mathbf{z}}_{{\tau }_{k+1}})$$

$$\forall B\in {\mathcal{E}}_{{\mathbf{\tau }}_{\mathbf{k}+1}},{\eta }_{k+1}(B)$$

$$Q_k(f)({\mathbf{Z}}_{{\tau }_k})={\int }_{{\mathbf{E}}_{{\tau }_{k+1}}}f({\mathbf{z}}_{{\tau }_{k+1}}){\mathbf{V}}_k(d{\mathbf{z}}_{{\tau }_{k+1}}|{\mathbf{Z}}_{{\tau }_k})G_k({\mathbf{Z}}_{{\tau }_k})$$

$${\Psi }_k(\eta )(f)=\frac{\int G_k(\mathbf{z})f(\mathbf{z})d\eta (\mathbf{z})}{\eta (G_k)}.$$

$${\eta }_k\stackrel{{\Psi }_k}{\stackrel{\text{mbox}}{\to }}{\tilde{\eta }}_k\stackrel{.V_k}{\stackrel{\text{mbox}}{\to }}{\eta }_{k+1}.$$

Full text

Used on such concentrated PDMPs, the IPS is inefficient and does not always provide a variance reduction. Indeed, the efficiency of the IPS method relies on simulating many different trajectories during its propagation steps, but unfortunately concentrated PDMPs are likely to generate the same deterministic trajectories over and over. We propose an adaptation of the IPS method called IPS+M that reduces this phenomenon. The IPS+M consists in modifying the propagation steps of the IPS, by conditioning the propagation to avoid generating the same trajectories multiple times. We prove that, compared to the IPS, the IPS+M method always provides an estimator with a lower variance. We also carry out a quick simulation study on a two-components system that confirms this result.

††journal: Chaos: An Interdisciplinary Journal of Nonlinear Science

1 Introduction

For both safety and regulation issues, the reliability of industrial systems has to be assessed. The considered systems (dams or nuclear power plants, for instance) are complex dynamic hybrid systems, so only simulation methods can be reasonably considered to assess their reliability. The failure of such a dynamic hybrid system generally corresponds to a physical variable of the system (temperature, pressure, water level) entering a critical region. The simulation of such a system requires to accurately model the trajectory of the physical variables. The evolution of these physical variables are generally determined by simple differential equations derived from the laws of physics. As the physical context giving the differential equations generally depends on the statuses of the multiple components of the systems (on, off, or failed), the differential equations can change whenever there is a change in the components statuses. To encounter for this hybrid interplay between the discrete process of components’ statuses and the continuous evolution of the physical variables, we model the evolution of the state of a system by a piecewise deterministic Markov process (PDMP) [5, 6, 7, 3]. PDMPs are meant to represent a large class of Markovian processes that do not include diffusion, as such, they benefit from high modeling capacity: they can model many industrial systems. For instance, EDF has recently developed the PyCATSHOO toolbox [2], which allows the modeling of dynamic hybrid systems, the main option within this toolbox is to evaluate the dependability criteria of the studied systems by Monte Carlo simulations.

As the industrial systems are highly reliable, their failure is a rare event, and the naive Monte-Carlo simulation method (MC) is computationally too intensive in this context. The objective of our work is to set up new algorithms to accelerate the reliability assessment of such industrial systems. To do so we want to use faster methods, such as variance reduction methods. A variance reduction method is a method that yields an unbiased estimator with a smaller variance than the Monte-Carlo estimator. The estimation being more accurate, we need less simulation runs to reach the desired accuracy, thus we save computational time. The variance reduction is generally achieved by altering the simulation process, and then correcting the bias induced by this modification of simulation process by appropriately weighting each simulation. Several acceleration methods for variance reduction can be proposed, such as importance sampling methods or particle filter methods (also called subset simulations). In this article we focus on one example of particle filter method: the interacting particle system (IPS) method [9].

Unlike in importance sampling methods, in the IPS method we keep simulating according to the original system. Another difference is that we do not simulate directly the trajectories on the entire observation duration, the observation duration is subdivided into small intervals of time and we simulate the trajectories sequentially one interval of time after the other one. These sequential simulations of the trajectories consist alternating between an exploration step and a selection step. During the exploration step we simulate the trajectories on a small time interval, therefore exploring the most probable trajectories on a short horizon of time. Then we apply a selection step on these trajectories replicating the trajectories which seem ”close” to failure and giving up the less ”promising” trajectories. At the next exploration step, only replicated trajectories are continued, before the next selection, and the next exploration and so on… until each selected trajectory reaches the end of the observation time window. This way the effort of simulation is concentrated on selected trajectories, which have higher chance of becoming a failing trajectory before the end of the observation time window. In the end we get more failing trajectories to fuel our estimation, and, if the selection was well done, the IPS method yields an unbiased estimator with a smaller variance than the MC estimator.

When we try to apply the IPS method in order to estimate the failure probability of such reliable complex hybrid systems, the IPS methods turn out to be inefficient. Indeed the estimation provided by the IPS method often has, in this case, a higher variance than the one of the Monte-Carlo estimator. IPS does not perform well in this context because the application case (a reliable complex hybrid system) makes it hard to conduct the exploration steps of the IPS method efficiently. Indeed, an industrial system is often modeled by what we call a ”concentrated PDMP”, which is a PDMP with low jump rates and concentrated jump kernels on boundaries. The typical jump rate is low because it is the sum of the failure rates of the components in working condition and of the repair rates of the failed components which are all very low. These failure and repair rates are very low because the components are reliable and their repairs are slow. The typical jump kernel is concentrated because most of the probability mass of a jump is concentrated on one (safe) output. Indeed the jump kernel on boundaries model the automatic control mechanisms within the system. During such a control mechanism there is a small probability that some component(s) fail on demand but the most likely output is that the system jumps on the safe state aimed by the control mechanism, consequently the probability mass of the jump kernel is concentrated on this output. Due to these characteristics of the model there is a high probability that no component failure or repair occurs during the short exploration time, and with a PDMP it means that all the trajectories are likely to follow the same deterministic paths. So when we explore the trajectory space most simulated trajectories end up being the same one, hence limiting our exploration of the trajectory space. To avoid this pitfall, a particular filter was proposed in [20] that enhances the occurrence of random jumps (failure or repairs) or modifies the occurrence time of the last jump. However the proposed method is limited to a different case of PDMP. This class of PDMP does not include concentrated PDMP, as it contains only PDMP without boundaries which allows continuous jumps kernels, i.e. it does not allow to model automatic control mechanisms in components. Moreover, an a priori bound on the number of jumps in a time interval is required, but in our case we do not have such information.

In order to adapt the IPS to concentrated PDMPS, we propose instead to use an approach based on the memorization method developed in [16]. The idea is to start the exploration by finding the most likely trajectories continuing each batch of replicated trajectories. Then, we condition the rest of the exploration to avoid these trajectories. As a result, the simulated trajectories have much more chances to differ which improves the quality of the exploration and reduce the variance of the estimator. To correct the bias induced by this modification of the simulation process, we have to modify the weight of each trajectory. We call our adaptation of the IPS to PDMP the IPS+M for sequential Monte-Carlo sampler with memorization method.

The rest of the paper paper is organized as follows: Section 2 is dedicated to the presentation of our model of the system, Section 3 presents the IPS method and introduces the optimal potential functions for this method, Section 5 presents the IPS+M method which adapts the IPS algorithm for PDMPs with low jump rates, Section 6 explains how to force the differentiation of the trajectories using the memorization method, and finally in Section 7 we illustrate the better efficiency of the IPS+M method on a toy example.

2 A model of the system based on a PDMP with discrete jump kernel

2.1 The model

We denote by $Z_{t}$ the state of the system at time $t$. $Z_{t}$ is the combination of the physical variables of the system, noted $X_{t}$, and of the statuses of all the components within the system, noted $M_{t}$: $Z_{t}=(X_{t},M_{t})$. We consider that $X_{t}\in\mathbb{R}^{d}$, and that $M_{t}\in\mathbb{M}=\{On,Off,F\}^{N_{c}}$ where $F$ corresponds to a failed status, and $N_{c}$ is the number of components in the system. The value of $M_{t}$ is sometimes referred as the mode of the system. Here we consider only three categories for the status of a component : $On,Off,F$, but it is possible to include more categories as long as the set of the possible modes $\mathbb{M}$ stays countable.

The process $Z_{t}$ is piecewise continuous, and each discontinuity is called a jump. Between two jumps there is no change in the components’ statuses, and the dynamics of the physical variables can be expressed thanks to an ordinary differential equation derived from the law of physics:

[TABLE]

We note $\phi_{(x,m)}(t)$ the solution of this equation when $X_{0}=x$ and $M_{0}=m$. Then for any time $s>0,\,$ if $T$ is the time separating $s$ from the next jump time, we have

[TABLE]

Similarly, a flow function on the states can be defined. If $z=(x,m)$, then we define \Phi_{z}(t)=\big{(}\phi_{(x,m)}(t),m\big{)}, and so

[TABLE]

As the physical variables are often continuous, the jumps are essentially used to model changes in the statuses of the components. These jumps can occur for two reasons.

Firstly, a jump can correspond to an automatic control mechanism (See Figure 1). In a given mode $m$, such mechanism is typically triggered when the physical variables cross some threshold. For each mode $m$, we define an open and connected set $\Omega_{m}$, so that these thresholds determine its boundary $\partial\Omega_{m}$. So when $M_{t}=m$ the value of the physical variables $X_{t}$ is restricted to the set $\Omega_{m}\subset\mathbb{R}^{d}$. Letting $E_{m}=\{(x,m),x\in\Omega_{m}\}$ be the set of the possible states with mode $m$, in terms of state a jump associated to a control mechanism is triggered whenever the state $Z_{s+t}$ hits the boundary of $E_{m}$. The set of possible states is therefore defined by:

[TABLE]

Secondly, jumps can correspond to a spontaneous failure or repair. In such case the jump occurs before $Z_{s+t}$ hits the boundary of $E_{m}$.

The occurrence time of a spontaneous jump is modeled by a jump rate $\lambda(Z_{s+t})$. For instance, when we consider that components fail or are repaired one at a time, this jump rate is the sum of the failure rates of the components in working condition and of the repair rates of the broken components. In a more general case, each transition from a mode $m$ to a mode $m^{+}$ has its own rate $\lambda_{m\to m^{+}}(X_{s+t})$, and we have

[TABLE]

We define the cumulative jump rate $\Lambda_{z}(t)$ by \Lambda_{z}(t)=\int_{0}^{t}\lambda\big{(}\Phi_{z}(u)\big{)}du, so that \forall t\in[0,T),\ \Lambda_{Z_{s}}(t)=\int_{0}^{t}\lambda\big{(}Z_{s+u}\big{)}du. Eventually the cumulative distribution function (cdf) of $T$ (the time until the next jump starting form a state $Z_{s}=z$) takes the form:

[TABLE]

Here $t^{*}_{z}=\inf\{t>0,\Phi_{z}(t)\in\partial E_{m}\}$ is the time until the flow hits the boundary starting from a state $z=(x,m)$. When there is no boundary i.e. $\{s>0,\Phi_{z}(s)\notin E_{m}\}=\emptyset$, we take the convention that t^{*}_{z}=+\infty\. With this definition, the law of $T$ has a continuous part associated to spontaneous jumps (spontaneous failures and repairs) and a discrete part associated to forced jumps (control mechanisms). A possible reference measure for $T$ knowing $Z_{s}=z$ is then

[TABLE]

where $leb(.)$ corresponds to the Lebesgue measure.

If a jump occurs at time $S$, then the distribution of the destination of the jump is expressed by a transition kernel $\mathcal{K}_{Z_{S}^{\text{-}}}$ where $Z_{S}^{\text{-}}$ is the departure state of the jump. If  \overline{\hbox{}}$$E\, is the closure of $E$, we have Z_{S}^{\text{-}}\in\hbox to0.0pt{\hskip 2.31245pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorE,}}. Let $Z_{S}^{+}\in E$ be the arrival state, and $\mathscr{B}(E)$ be the Borelian $\sigma$-algebra on $E$. Then, the law of a jump from a departure state $z^{\text{-}}$ is defined by:

[TABLE]

For any departure state z^{\text{-}}\in\hbox to0.0pt{\hskip 2.31245pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorE,}}, we define $x_{a}(z^{\text{-}},m^{+})$ as the arrival for the vector of physical variables when the transition $m^{\text{-}}\to m^{+}$ is triggered from a state $z^{\text{-}}$. We define by z_{a}(z^{\text{-}},m^{+})=\big{(}x_{a}(z^{\text{-}},m^{+}),m^{+}\big{)} the arrival state when the transition $m^{\text{-}}\to m^{+}$ is triggered from a state $z^{\text{-}}$. Using this definition, for any departure state z^{\text{-}}\in\hbox to0.0pt{\hskip 2.31245pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorE,}}, we define the measure $\nu_{z^{\text{-}}}$ on \big{(}E,\mathscr{B}(E)\big{)} by:

[TABLE]

The kernel $\mathcal{K}_{z^{\text{-}}}$ is absolutely continuous with respect to the measure $\nu_{z^{\text{-}}}$. Denoting by $K_{z^{\text{-}}}$ its density with respect to $\nu_{z^{\text{-}}}$ we have:

[TABLE]

Note that we only consider systems for which $\nu_{z^{\text{-}}}$ is discrete, whatever the departure state $z^{\text{-}}\in E$ may be. This hypothesis of discreteness of $\nu_{z^{\text{-}}}$ is mandatory to apply the method presented in Section 5.

If $\nu_{z}$ has a Dirac contribution at point $z$, then the kernel must satisfy $K_{z}(z)=0$ , so it is not possible to jump on the departure state. In some applications the physical variables are continuous, so that $K_{z^{\text{-}}}(z)$ is zero whenever $x^{-}\neq x$. In some cases, one might want the model to include renewable and aging components. Then the vector $X_{t}$ should include the time since the last renewal of such component, and so the vector $X_{t}$ can be discontinuous at the time of a renewal. In both situations, when $z^{\text{-}}$ is not on the boundary of $E_{m^{\text{-}}}$, the jump kernel has this form:

[TABLE]

When $z^{\text{-}}\in\partial E_{m^{\text{-}}}$, a control mechanism is triggered : some components are required to turn on, or to turn off, so that the system reaches a desired state $z_{c}$. The transition to this state is usually very likely, so when $z^{\text{-}}\in\partial E_{m^{\text{-}}}$, the jump kernel $K_{z^{\text{-}}}$ tends to concentrate all or a big portion of the probability mass on this state $z_{c}$. In industrial systems there can be components which have a small probability to fail when they are asked to turn on. This phenomenon is referred as a failure on demand. We denote by $\kappa_{i}(z^{-})$ the probability that the $i^{th}$ component fails on demand when the control is triggered. We denote by $fod(i,m^{\text{-}},m)$ the indicator being equal to one if the $i^{th}$ component fails on demand during the transition from $m^{\text{-}}$ to $m$, and to zero otherwise. We also define $ask(z^{\text{-}},z^{+})$ as the set gathering the indices of the components supposed to turn on during the control mechanism triggered from $z^{\text{-}}$ to $z^{+}$. Finally, when $z^{\text{-}}\in\partial E_{m^{\text{-}}}$, we have that:

[TABLE]

To generate $\mathbf{Z}_{t}=(Z_{s})_{s\in[0,t]}$ a trajectory of the states of the system, one can repeat the following steps: Starting with $s=0$,

Given a starting state $Z_{s}=z$, generate $T$ the time until the next jump using (3), 2. 2.

Follow the flow $\Phi$ until $s+T$ using (1), and set the departure state of the next jump as being $Z_{T}^{\text{-}}=\Phi_{z}(T)$, 3. 3.

generate $Z_{s+T}$ the arrival state of the next jump using $K_{Z_{T}^{\text{-}}}$ 4. 4.

repeat starting with $s=s+T$ until one gets a trajectory of size $t$.

Defined in this way, the process $Z_{t}$ is Markovian [4].

2.2 Law of trajectory

Let us denote by $n(\mathbf{Z}_{t})$ the number of jumps in the trajectory $\mathbf{Z}_{t}$, and by $S_{k}$ the time of the $k^{th}$ jump in the trajectory $\mathbf{Z}_{t}$ (with the convention that $S_{0}=0$ and $S_{n(\mathbf{Z}_{t})+1}=t$). $T_{k}=S_{k+1}-S_{k}$ denotes the time between two consecutive jumps. We define the $\sigma$-algebra $\mathscr{S}_{t}$ as the $\sigma$-algebra generated by the sets in \underset{\ n\in\mathbb{N}^{*}}{\bigcup}\mathscr{B}\Big{(}\Big{\{}\big{(}z_{k},t_{k}\big{)}_{k\leq n}\in(E\times\mathbb{R}^{*}_{\text{+}})^{n},\,\overset{n}{\underset{i=0}{\sum}}t_{i}=t\Big{\}}\Big{)}, where $\mathscr{B}(.)$ indicates the Borelians of a set. Letting $\Theta_{t}:\mathbf{Z}_{t}\to\left((Z_{k},T_{k})\right)_{0\leq k\leq n(\mathbf{Z}_{t})}$ be the application giving the skeleton of a trajectory, we can define the law of trajectories as an image law through $\Theta_{t}$. We have that for $B\in\mathscr{S}_{t}$:

[TABLE]

where $z_{j}^{-}=\Phi_{z_{j-1}}(t_{j-1})$, and $t^{*}_{n}=t-\sum_{i=0}^{n-1}t_{i}$.

In the rest of the paper we will denote $\mathbf{E}_{t}$ the set of trajectories of size $t$ that satisfy (1), and by $\mathscr{M}(\mathbf{E}_{t})$ the set of bounded $(\Theta_{t}^{-1}(\mathscr{S}_{t}),\mathbf{E}_{t})$-measurable functions.

2.3 A hybrid reference measure for the trajectories of PDMP

As they involve hybrid reference measures for the times between jumps, PDMPs are very degenerate processes. As a result, some of their realizations have a strictly positive probability to be simulated, but some other do not. The trajectories that concentrate a part of the probability mass on their own, i.e. that verify \mathbb{P}\big{(}\mathbf{Z}_{t}=\mathbf{z}_{t}\big{)}>0, are called the preponderant trajectories. For example, assuming that $\lambda$ is positive, if the trajectory $\mathbf{z}_{t}$ involves no jump, it is preponderant because we have \mathbb{P}\big{(}\mathbf{Z}_{t}=\mathbf{z}_{t}\big{)}=1-\exp\big{(}-\Lambda_{z_{0}}(t)\big{)}>0. Also, as we consider that the reference measures $\nu_{z^{\text{-}}}$ are discrete, there can be an other type of preponderant trajectories, which are the trajectories that only jump on the boundaries $\partial E_{m}$. Indeed, such trajectories $\mathbf{z}_{t}$ would verify

[TABLE]

where $s_{k}$ is the time of the $k^{th}$ jump in $\mathbf{z}_{t}$. Conversely, some realizations can be considered as negligible as they verify \mathbb{P}\big{(}\mathbf{Z}_{t}=\mathbf{z}_{t}\big{)}=0. These trajectories are the ones that involve a spontaneous jump, i.e. a jump starting from the interior of a set $E_{m}$.

This can be better understood by looking at the equation (11) which shows that the measure $\zeta_{z_{0},t}$ defined by :

[TABLE]

where $z_{j}^{-}=\Phi_{z_{j-1}}(t_{j-1})$, and $t^{*}_{n}=t-\sum_{i=0}^{n-1}t_{i}$ is a reference measure for the law of trajectories. This reference measure highlights the fact that the trajectories with no jump or with only jumps on boundaries can concentrate some probability mass on their own, because they refer to a Dirac contribution of the measure $\zeta_{z_{0},t}$. Indeed for such trajectories the times between two consecutive jumps verify $t_{k}=t_{z_{s_{k}}}^{*}$ $\forall k<n$, and so, they always refer to the discrete part of the measures $\mu_{z_{k-1}}$ $\forall k<n$, therefore such trajectories are related to the discrete part of $\zeta_{z_{0},t}$. Equation (13) also shows that the remaining probability mass is distributed continuously among the trajectories with at least one spontaneous jump. Indeed, in a negligible trajectory, if for example the $k+1^{th}$ jump is a spontaneous jump, then $t_{k}$ relates to the continuous part of the reference measure $\mu_{z_{k-1}}$.

2.4 Concentrated PDMP

For reliability assessment of a highly reliable system one often models the system by a PDMP with low jump rates and concentrated jump kernels on the boundaries. Indeed, the components of the system are often reliable and their repair takes time, hence the low jump rates, and as failures on demand during a control mechanism are unlikely the jump kernels on boundaries are concentrated on one safe arrival state (i.e. the state aimed by the control mechanism). We call this kind of PDMP a concentrated PDMP, mainly because the law of one trajectory concentrates a big part of its probability mass.

This trajectory happens to be the trajectory with no failure and no repair. As jump rates are low the probability of not having a spontaneous jump is close to one. For instance at a $k+1$-th jump this probability verifies

[TABLE]

So only jumps on boundaries are likely, and when the process hits a boundary $\partial E_{m}$, the arrival state aimed by the control mechanism is very likely. Denoting by $z_{s_{k}}$ this state for a $k$-th jump we have :

[TABLE]

So if $\mathbf{z}_{t}$ is a trajectory with no failure and no repair of reasonable size we have:

[TABLE]

In this article we adapt the IPS method for concentrated PDMPs. It is interesting to make the connection between our work and the modified particle filter developed in [12] which can be applied to the related context of discrete time processes with discrete spaces.

2.5 Reliability assessment

Let $t_{f}$ be an observation time. We denote by $\mathbf{E}_{t_{f}}$ the set of trajectories of size $t_{f}$ that verify (1), and $h\in\mathscr{M}(\mathbf{E}_{t_{f}})$. The methods presented in this article can be used for the estimation of any quantity $p_{h}$ defined by:

[TABLE]

We are interested in estimating the probability, noted $p_{\mathscr{D}}$, that the system fails before the final observation time $t_{f}$ knowing it was initiated in a state $z_{0}$. Letting $\mathscr{D}$ be the set of trajectories of length $t_{f}$ which pass through the critical region $D\subset E$, we have

[TABLE]

Although our application relates to the case where $h=\mathbbm{1}_{\mathscr{D}}$, the IPS method will be presented with an arbitrary (bounded) function $h$.

3 The IPS method

3.1 A Feynman-Kac model

For a measure $\eta$ and a bounded measurable function $f$ we note $\eta(f)=\int f\,d\eta$. For a measure $\eta$ and a kernel $V$, $\eta V$ denotes the measure such that $\eta V(h)=\int\int h(y)V(dy|x)d\eta(x)$, and for a bounded measurable function $f$, $V(f)$ is the function such that $V(f)(x)=\int h(y)V(dy|x)$.

Consider a subdivision of the interval $[0,t_{f}]$ into $n$ sub-intervals of equal lengths, noted $[\tau_{k},\tau_{k+1})$, and such that $0=\tau_{0}<\tau_{1}<\dots<\tau_{n\,\text{-}1}<\tau_{n}=t_{f}$. Let $V_{k}$ be the Markovian transition measure extending a trajectory of size $\tau_{k}$ into a trajectory of size $\tau_{k+1}$, such that

[TABLE]

For each $k<n$ we denote $G_{k}$ the potential function on $\mathbf{E}_{\tau_{k}}$, such that:

[TABLE]

These potential functions are the main inputs of the method. The choice of the potential functions is important, because it will ultimately determine the variance of the estimator of $p_{h}$ provided by the IPS method. A good choice depends on the system and on the target function $h$. The potential functions are used to define the target probability measures $\tilde{\eta}_{k}$ for each $k\leq n$ , such that:

[TABLE]

or equivalently

[TABLE]

Originally the IPS method comes from filtering methods. Filtering methods aim at estimating the target measures and, in these methods, the potential functions are chosen so that the $\eta_{k}$ match the target measures. But in our cases we have no interest in estimating the target measures, we only want to estimate $p_{h}$. So the potential functions can be chosen more freely. They are used to propose a probabilistic representation of $p_{h}$ in terms of a selection+mutation dynamics, which makes it possible to build an estimator of $p_{h}$ with a reduced variance.

We define the propagated target measures $\eta_{k}$ such that $\eta_{0}=\tilde{\eta}_{0}$ and for $k\geq 0$, $\eta_{k+1}=\tilde{\eta}_{k}V_{k}$. We have :

[TABLE]

or equivalently

[TABLE]

For $k=0$ we consider that $\eta_{0}=\delta_{0}$, but the methods would still be valid if we had $\eta_{0}\neq\delta_{0}$. We define $Q_{k}$ such that for $f\in\mathscr{M}(\mathbf{E}_{\tau_{k+1}})$,

[TABLE]

and set $Q_{k,n}=Q_{k}Q_{k+1}\dots Q_{n}$. Let $\Psi_{k}$ be the application that transforms a measure $\eta$ defined on $\mathbf{E}_{\tau_{k}}$ into a measure $\Psi_{k}(\eta)$ defined on $\mathbf{E}_{\tau_{k}}$ as follows:

[TABLE]

We say that $\Psi_{k}(\eta)$ gives the selection of $\eta$ through the potential $G_{k}$. Notice that $\tilde{\eta}_{k}$ is the selection of $\eta_{k}$ as $\tilde{\eta}_{k}=\Psi_{k}(\eta_{k})$. The target distributions can therefore be built according to the following pattern, in which a propagation step follows a selection step:

[TABLE]

We also define the associated unnormalized measures $\tilde{\gamma}_{k}$ and $\gamma_{k+1}$, such that for $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$:

[TABLE]

and for $f\in\mathscr{M}(\mathbf{E}_{\tau_{k+1}})$:

[TABLE]

Denoting $f_{h}(\mathbf{Z}_{\tau_{n}})=\frac{h(\mathbf{Z}_{\tau_{n}})}{\prod_{s=0}^{n-1}G_{s}(\mathbf{Z}_{\tau_{s}})}$, notice that we have:

[TABLE]

3.2 The IPS algorithm and its estimators

The IPS method provides an algorithm to generate weighted samples which approximate the probability measures $\eta_{k}$ and $\tilde{\eta}_{k}$ respectively for each step $k$. For the sample approximating $\eta_{k}$, we denote $\mathbf{Z}_{\tau_{k}}^{j}$ the $j^{th}$ trajectory and $W_{k}^{j}$ its weight. Respectively, for the sample approximating $\tilde{\eta}_{k}$, we denote $\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{j}$ the $j^{th}$ trajectory and $\accentset{\sim}{W}_{k}^{j}$ its associated weight. For simplicity reasons, in this paper, we consider that the samples all contain $N$ trajectories, but it is possible to modify the sample size at each step, as illustrated in [17]. The empirical approximations of $\eta_{k}$ and $\tilde{\eta}_{k}$ are denoted by $\eta_{k}^{N}$ and $\tilde{\eta}_{k}^{N}$ and are defined by:

[TABLE]

So for all $k\leq n$ and $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$,

[TABLE]

By plugging these estimations into equations (22) and (22), we get estimations for the unnormalized distributions. Denoting by $\tilde{\gamma}_{k}^{N}$ and $\gamma_{k}^{N}$ these estimations, for all $k\leq n$ and $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$, we have:

[TABLE]

Plugging the estimations $\eta_{k}^{N}$ into equation (24), we get an estimator $\hat{p}_{h}$ of $p_{h}$ defined by:

[TABLE]

The algorithm builds the samples sequentially, alternating between a selection step and a propagation step. The $k^{th}$ selection step transforms the sample $(\mathbf{Z}_{k}^{j},W_{k}^{j})_{j\leq N}$ into the sample $(\accentset{\sim}{\mathbf{Z}}_{k}^{j},\accentset{\sim}{W}_{k}^{j})_{j\leq N}$. This transformation is done with a multinomial resampling scheme: the $\accentset{\sim}{\mathbf{Z}}_{k}^{j}$’s are drawn with replacement from the sample $(\mathbf{Z}_{k}^{j})_{j\leq N}$, each trajectory $\mathbf{Z}_{k}^{j}$ having a probability $\frac{W_{k}^{j}G_{k}(\mathbf{Z}_{\tau_{k}}^{j})}{\sum_{i=1}^{N}W_{k}^{i}G_{k}(\mathbf{Z}_{\tau_{k}}^{i})}$ to be drawn each time. We denote by $A^{j}_{k}$ the ancestor index of the $j^{th}$ trajectory in the selected sample, such that $\accentset{\sim}{\mathbf{Z}}_{k}^{j}=\mathbf{Z}_{\tau_{k}}^{A^{j}_{k}}$. We let $\tilde{N}_{k}^{j}=card\{i,A^{i}_{k}=j\}$ be the number of times the particle $\mathbf{Z}_{k}^{j}$ is replicated in the sample $(\accentset{\sim}{\mathbf{Z}}_{k}^{j},\accentset{\sim}{W}_{k}^{j})_{j}$, so $N=\sum_{j=1}^{N}\tilde{N}_{k}^{j}$. After this resampling the weights $\accentset{\sim}{W}_{k}^{j}$ are set to $\frac{1}{N}$. The interest of this selection by resampling is that it discards low potential trajectories and replicates high potential trajectories. So the selected sample focuses on trajectories that will have a greater impact on the estimations of the next distributions once extended.

If one specifies potential functions that are not positive, there can be a possibility that at a step $k$ we get $\forall j,G_{k}(\mathbf{Z}_{\tau_{k}}^{j})=0,$ and so the probability for resampling cannot be defined. When this is the case, the algorithm stops and we consider that $\forall s\geq k$ the measures $\tilde{\eta}_{s}^{N}$ and $\eta_{s}^{N}$ are equal to the null measure.

Then the $k^{th}$ propagation step transforms the sample $(\accentset{\sim}{\mathbf{Z}}_{k}^{j},\accentset{\sim}{W}_{k}^{j})_{j\leq N}$, into the sample $(\mathbf{Z}_{k+1}^{j},W_{k+1}^{j})_{j\leq N}$. Each trajectory $\mathbf{Z}_{k+1}^{j}$ is obtained by extending the trajectory $\accentset{\sim}{\mathbf{Z}}_{k}^{j}$ on the interval $[\tau_{k},\tau_{k+1})$ using the transition kernel $V_{k}$. The weights satisfy $W_{k+1}^{j}=\accentset{\sim}{W}_{k}^{j},\ \forall j$. Then the procedure is iterated until the step $n$. The full algorithm to build the samples is displayed in Figure 3.

For the sake of simplicity, we will make the following assumption: $\exists\,\varepsilon_{1},\varepsilon_{2}\in\mathbb{R}^{+}$ such that $\forall\,\mathbf{Z}_{\tau_{k}}\in\mathbf{E}_{\tau_{k}}$:

[TABLE]

Theorem 1**.**

When (G) is verified the estimator (28) is unbiased and strongly consistent.

The proof of theorem 1 follows from Theorems 7.4.2 and 7.4.3 in [8].

Theorem 2**.**

When (G) is verified:

[TABLE]

*where *

[TABLE]

A proof of this CLT can be found in [8] chapter 9 at the theorem 9.3.1 . For the estimation of the variance $\sigma^{2}_{IPS,G}$ we refer the reader to [17].

3.3 Classical improvements of the IPS method: the SMC method, and its alternative resampling

We have seen that the resampling steps have the advantage of replicating high potential trajectories and discarding low potential trajectories. However the resampling steps also introduce some additional fluctuations into the estimation (see (31)). So we would like to trigger them only when it is judicious. Typically, not when the potentials of all the trajectories are similar, as in this case there is not point in discarding or replicating some trajectories over others. In order to avoid pointless resampling, one can trigger the selection step only when the weights are unbalanced. This is done in the Sequential Monte Carlo (SMC) algorithm with adaptive resampling presented in Figure 4. In this algorithm, the heterogeneity of the weights is quantified using the effective sample size. At the $k^{th}$ step the effective sample size is defined by:

[TABLE]

Its value is between $1$ and $N$ and it measures the homogeneity in the candidate weights $\frac{W_{k}^{i}G_{k}(\mathbf{Z}_{\tau_{k}}^{i})}{\sum_{j}W_{k}^{j}G_{k}(\mathbf{Z}_{\tau_{k}}^{j})}$: when $ESS_{k}=N$ the weights are perfectly balanced and are all equal to $\frac{1}{N}$, and conversely when $ESS_{k}=1$ all the weights are null except one, which concentrates the totality of the mass. Therefore, one considers the weights are too unbalanced when $ESS_{k}\leq eN$ where $e\in[0,1]$ is a tuning parameter.

Note that in the presented algorithm one can use alternative strategies to select high potential trajectories. Here, the presented algorithms include a standard multinomial resampling procedure, but one can also use residual resampling or stratified resampling without altering the properties of the estimator. Empirical results suggest that these alternative resampling schemes yield estimations with smaller variances [11, 14]. There are also recent theoretical results on the performance of different resampling schemes [13]. MCMC steps with invariant distribution $\tilde{\eta}_{k}$ can also be included in the algorithm after the resampling step. Some adaptations of the algorithm for parallel implementations have also been studied in [19] for instance.

4 Choice of the potential functions

Note that the variance of $\hat{p}_{h}$ depend on the number of subdivisions and on the choice of the potential functions. We display here an important hint on how to select potential functions that yield a small variance. Indeed, the theoretical expressions of the potential functions that minimize the asymptotic variance of the IPS estimator are known [1].

Theorem 3**.**

For $k\geq 1$, let $G^{*}_{k}$ be defined by:

[TABLE]

if \mathbb{E}\Big{[}\mathbb{E}\big{[}h(\mathbf{Z}_{\tau_{n}})\big{|}\mathbf{Z}_{\tau_{k}}\big{]}^{2}\big{|}\mathbf{Z}_{\tau_{k-1}}=\mathbf{z}_{\tau_{k-1}}\Big{]}\neq 0, and $G^{*}_{k}(\mathbf{z}_{\tau_{k}})=0$ otherwise. For $k=0$, we define

[TABLE]

The potential functions minimizing $\sigma^{2}_{IPS,G}$ are the ones that are proportional to the $G^{*}_{k}$’s $\forall k\leq n$. The optimal variance of the IPS method with $n$ steps is then

[TABLE]

For the reliability assessment case, the optimal target distributions have the form:

[TABLE]

and the potential functions would be defined by (33) and (34), taking $h=\mathbbm{1}_{\mathscr{D}}$. As we do not have closed-form expressions of the functions \mathbf{z}_{\tau_{k}}\to\sqrt{\mathbb{E}\big{[}\mathbb{E}[\mathbbm{1}_{\mathscr{D}}(\mathbf{Z}_{\tau_{n}})|\mathbf{Z}_{\tau_{s+1}}]^{2}|\big{|}\mathbf{Z}_{\tau_{s}}=\mathbf{z}_{\tau_{k}}\big{]}}, we propose to use instead a parametric approximation of these expectations based on our knowledge of the system and denoted by $U_{\alpha}(\mathbf{Z}_{\tau_{s}})$ so that we take

[TABLE]

For a system including similar components in parallel redundancy, we propose to set $U_{\alpha}(\mathbf{Z}_{\tau_{s}})=1$ when $\mathbf{Z}_{\tau_{s}}$ has already reached the failure region once, and to set U_{\alpha}(\mathbf{Z}_{\tau_{s}})=\exp\big{[}-\alpha\,(b(Z_{\tau_{s}})+1)^{2}\big{]}\,L(\tau_{s}) otherwise, where $L$ is a positive function, and $b(Z)$ indicates the number of working components within a state $Z$. Here $\alpha$ is a parameter tuning the strength of the selection.

5 The IPS+M method for concentrated PDMPs

5.1 The issue with PDMP modeling a reliable system

When it is used on a reliable system and therefore on a concentrated PDMP (see Section 2.4), the IPS method tends to loose in efficiency. This efficiency loss can be attributed to the exploration steps. Remember that an exploration step comes after a selection step: it builds a sample $(\mathbf{Z}_{\tau_{k+1}}^{j},W_{k+1}^{j})_{j\leq N}$ by extending the trajectories of a selected sample $(\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{j},\accentset{\sim}{W}_{k}^{j})_{j\leq\tilde{N}}$. This newly built sample $(\mathbf{Z}_{\tau_{k+1}}^{j},W_{k+1}^{j})_{j\leq N}$ fulfills two goals : 1) It contributes to the empirical approximation $\eta_{k+1}^{N}$ of $\eta_{k+1}$. 2) It is used as a candidate sample for the next selection. But this second goal is often poorly achieved with a concentrated PDMP. Indeed, in order to get a good approximation $\tilde{\eta}_{k+1}^{N}$ of $\tilde{\eta}_{k+1}$, it is preferable that the candidate sample to selection $(\mathbf{Z}_{\tau_{k+1}}^{j},W_{k+1}^{j})_{j\leq N}$ contains as many different trajectories as possible, along with high potential trajectories. Unfortunately, with this kind of PDMP, it is generally not the case: the candidate sample often contains several replicates of the same trajectories, and no high potential trajectory. Therefore each distribution $\tilde{\eta}_{k}$ is poorly represented, and so is each target distribution $\eta_{k+1}$, which eventually deteriorates the quality of the estimator $\hat{p}_{\mathscr{D}}$ or $\hat{p}_{h}$.

To understand why the exploration steps are not likely to generate many different trajectories with a concentrated PDMP, we have to come back at the beginning of the propagation step. At that point, the sample $(\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{j},\accentset{\sim}{W}_{k}^{j})_{1\leq j\leq N}$ is naturally clustered because of the previous selection step, each of the clusters containing several replicates of the same trajectories. We can rewrite (25) in the following way

[TABLE]

where $\sum_{j=1}^{N}\tilde{N}_{k}^{j}=N$. In practice many of the $\tilde{N}_{k}^{j}$ are null and only a few are positive and the $N$ resampled trajectories $(\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{j})_{j\leq N}$ are concentrated on a few trajectories. Then, each of the $\tilde{N}_{k}^{j}$ trajectories of the $j$-th cluster is extended by using the same distribution $V_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$. (For all index $i$ such that $A_{k}^{i}=j$ the trajectory $\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{i}$ is extended with the kernel $V_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$). As the kernel $V_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$ corresponds to a concentrated PDMP, it is likely to extend all the trajectories of a cluster in the same manner. The trajectory $\mathbf{a}_{\tau_{k+1}}^{k,j}$ which extends $\mathbf{Z}_{\tau_{k}}^{j}$ until $\tau_{k+1}$ without spontaneous jump or failure concentrates the mass of the kernel $V_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$. Indeed, at this point we have :

[TABLE]

Therefore each of the trajectories $\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{i}=\mathbf{Z}_{\tau_{k}}^{j}$ in a cluster tends to be extended in $\mathbf{a}_{\tau_{k+1}}^{k,j}$. Thus, the trajectories within a cluster are likely to stay clumped together during the propagation, and the propagated sample $(\mathbf{Z}_{\tau_{k+1}}^{j},W_{k+1}^{j})_{j\leq N}$ is very likely to be clustered too. When the preponderant trajectories $\mathbf{a}_{\tau_{k+1}}^{k,j}$ have low potential values, the sample is not likely to contain high potential trajectories. Consequently the selection step having no good candidates and too few candidates, it tends to yield an inaccurate estimation of the distributions $\tilde{\eta}_{k}$.

This situation is typical of reliability assessment. In that context, a well constructed potential function is close to $G^{*}_{k}$ wherein $h=\mathbbm{1}_{\mathscr{D}}$. So the potential of a trajectory $G_{k+1}(\mathbf{Z}_{\tau_{k+1}})$ should be high if its final state $Z_{\tau_{k+1}}$ is more degraded than the state $Z_{\tau_{k}}$. This generally implies that $\mathbf{Z}_{\tau_{k+1}}$ includes at least one component failure between $\tau_{k}$ and $\tau_{k+1}$. As the preponderant trajectories $\mathbf{a}_{\tau_{k+1}}^{k,j}$ do not contain failure between $\tau_{k}$ and $\tau_{k+1}$ they generally are associated with low potential values.

The segment of $\mathbf{a}_{\tau_{k+1}}^{k,j}$ on $(\tau_{k},\tau_{k+1}]$ relates to a Dirac contribution of the measure $\zeta_{Z^{j}_{\tau_{k}},\tau_{k+1}-\tau_{k}}$. We can, therefore, decompose the expected propagation of the trajectory $\mathbf{Z}_{\tau_{k}}^{j}$ in this way:

[TABLE]

where $f\in\mathscr{M}(\mathbf{E}_{\tau_{k+1}})$. And the expected propagation of $\eta_{k}^{N}$ would be:

[TABLE]

5.2 Modify the propagation of clusters

In order to diversify the simulated trajectories, and to increase the precision of the estimation, we propose to modify the way we extend the selected trajectories. Here we consider that the size of the propagated sample can differ from the size of the previous selected sample. We now denote $\tilde{N}_{k}$ the size of the $k^{th}$ selected sample, and $N_{k+1}$ the size of the $k^{th}$ propagated sample, with the convention $N_{0}=N$. We stressed out, in section 5.1, that the propagation step aims at providing an estimation of $\eta_{k+1}=\tilde{\eta}_{k}V_{k}$ using the selected sample $(\accentset{\sim}{\mathbf{Z}}_{k}^{j},\accentset{\sim}{W}_{k}^{j})_{j\leq\tilde{N}_{k}}$. In other words, the selection step aims at providing a propagated weighted sample $(\mathbf{Z}_{\tau_{k+1}}^{j},W_{k+1}^{j})_{j\leq N_{k+1}}$ to estimate the distribution $\tilde{\eta}_{k}^{\tilde{N}_{k}}V_{k}$ defined by:

[TABLE]

where $f\in\mathscr{M}(\mathbf{E}_{\tau_{k+1}})$. We denote by $\bar{V}_{k}$ the Markovian kernel from $\mathbf{E}_{\tau_{k}}$ to $\mathbf{E}_{\tau_{k+1}}$ such that, for any trajectory $\mathbf{Z}_{\tau_{k}}^{j}$, $\bar{V}_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$ is the conditioning of $V_{k}(.|\mathbf{Z}_{\tau_{k}}^{j})$ to $\mathbf{E}_{k}\backslash\{\mathbf{a}_{\tau_{k+1}}^{k,j}\}$. $\bar{V}_{k}$’s density verifies:

[TABLE]

Using (43) we can decompose $\tilde{\eta}_{k}^{\tilde{N}_{k}}V_{k}$ as follows:

[TABLE]

In the original IPS algorithm, the sample approximating $\tilde{\eta}_{k}^{\tilde{N}_{k}}V_{k}$ is built by directly extending each trajectory in the selected sample. When we extend the replicates of a cluster, in average a proportion $V_{k}(\mathbf{a}_{\tau_{k+1}}^{k,j}|\mathbf{Z}_{\tau_{k}}^{j})$ of the replicates are extended in $\mathbf{a}_{\tau_{k+1}}^{k,j}$. This proportion of trajectories extended in $\mathbf{a}_{\tau_{k+1}}^{k,j}$ then serves as an estimation of $V_{k}(\mathbf{a}_{\tau_{k+1}}^{k,j}|\mathbf{Z}_{\tau_{k}}^{j})$. But it is not necessary to misspend all these replicates to estimate the probability of the preponderant trajectory. If we use equation (45), we would need to generate the trajectory $\mathbf{a}_{\tau_{k+1}}^{k,j}$ only once to assess its contribution to the propagation of the cluster. Also, $\mathbf{a}_{\tau_{k+1}}^{k,j}$ is easy to get. To generate it, it generally suffices to run the simulation process starting from the state $Z_{\tau_{k}}^{j}$ until time $\tau_{k+1}$, while setting the jumps rates and the probability of failure on demand to zero.

Therefore, for each cluster, we propose to use an additional replicate to generate $\mathbf{a}_{\tau_{k+1}}^{k,j}$ and compute exactly its contribution. So for any $j\in\{1,\dots N_{k}\}$, we will extend the selected trajectory $\mathbf{Z}_{k}^{j}$, $N_{k}^{j}$ times, where $N_{k}^{j}=\tilde{N}_{k}^{j}+\mathbbm{1}_{\tilde{N}_{k}^{j}>0}$. We denote $j_{i}$ the index of the $i^{th}$ replicates of $\mathbf{Z}_{k}^{j}$, and consider the added replicate has index [math] such that for $i\in\{0,\dots,\tilde{N}^{j}_{k}\}$ we have $\accentset{\sim}{\mathbf{Z}}^{j_{i}}_{\tau_{k}}=\mathbf{Z}_{k}^{j}$. The additional replicate is deterministically extended to the preponderant trajectory, so we have $\mathbf{Z}_{\tau_{k+1}}^{j_{0}}=\mathbf{a}_{\tau_{k+1}}^{k,j}$, and we set its weight to $W_{k+1}^{j_{0}}=V_{k}(\mathbf{a}_{\tau_{k+1}}^{k,j}|\mathbf{Z}_{\tau_{k}}^{j})\frac{\tilde{N}_{k}^{j}}{\tilde{N}_{k}}$, so that it carries all the mass associated to the preponderant trajectories of a cluster. Then we can use all the remaining $\tilde{N}_{k}^{j}$ trajectories in the cluster to estimate the non preponderant part of the cluster’s propagation (the $1^{st}$ term in the right hand side of equation (41)). For $i>0$, we condition the extensions to avoid $\mathbf{a}_{\tau_{k+1}}^{k,j}$ generating the $\mathbf{Z\,}_{\tau_{k+1}}^{j_{i}}$ according to the kernel $\bar{V}_{k}(.|\mathbf{Z\,}_{\tau_{k}}^{j})$ and set $W_{k+1}^{j_{i}}=\frac{1-V_{k}(\mathbf{a}_{\tau_{k+1}}^{k,j}|\mathbf{Z}_{\tau_{k}}^{j})}{\tilde{N}_{k}^{j}}\frac{\tilde{N}_{k}^{j}}{\tilde{N}_{k}}$.

Usually, the simulations of a restricted law are carried out using a rejection algorithm, but in our case a rejection algorithm would perform poorly. The rate of rejection would be too high, as it would be equal to $V_{k}(\mathbf{a}_{\tau_{k+1}}^{k,j}|\mathbf{Z}_{\tau_{k}}^{j})$ which is typically close to 1. For PDMPs, such simulations, conditioned to avoid a preponderant trajectory, can be efficiently carried out using the memorization method. This method, introduced in [16], shares similarities with the inverse method. It therefore benefits from not using any rejection, and so it is well suited to our applications. The memorization method is presented in section 6.

The target distributions $\tilde{\eta}_{k}$ are still estimated with $\tilde{\eta}_{k}^{\tilde{N}_{k}}$, using equation (39), but for $k=0$ to $n-1$, $\eta_{k+1}$, the propagation of a target distribution, is now estimated by :

[TABLE]

Let $\tilde{\mathbf{N}}_{k}=(N_{0},\tilde{N}_{0},N_{1},\tilde{N}_{1},\dots,\tilde{N}_{k})$ and $\mathbf{N}_{k}=(N_{0},\tilde{N}_{0},N_{1},\tilde{N}_{1},\dots,N_{k})$ We now note $\tilde{\gamma}_{k}^{\tilde{\mathbf{N}}_{k}}$ and $\gamma_{k}^{\mathbf{N}_{k}}$ the estimations of the unnormalized distributions, and for all $k\leq n$ and $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$, we define them by:

[TABLE]

In the end, $p_{h}$ is estimated using the equation :

[TABLE]

The full modified version of the algorithm is presented in Figure 5. We call this modified version of the IPS algorithm the IPS+M algorithm.

Throughout the rest of the paper, the notation $\mathbb{E}_{{}_{M}}$ will indicate that the expectation is associated to the IPS+M method and $\mathbb{E}$ will still denote the expectation for the original IPS method.

5.3 Convergence properties of the IPS+M estimators

In this section we show that the estimator $\hat{p}_{h}$ obtained with the IPS+M method have the same basic properties as the IPS estimator. With the IPS+M method, $\hat{p}_{h}$ converges almost surely to $p_{h}$, it is unbiased, and it satisfies a CLT. The proofs that we provide in this section follow the reasoning of the proofs in [8]. We present how to adjust the original proofs to take into account that the extensions of the trajectories within a cluster are no longer identically distributed. Finally we show that the asymptotic variance of the CLT is reduced with the IPS+M method.

5.3.1 The martingale decomposition of the anticipated biases

With the IPS+M method, we assumed that $\forall k,\,\tilde{N}_{k}=N$. For $p\leq 2n,$ we define $\mathcal{F}_{p}$ the filtration associated to the sequence of the $p$ first random samples built with the IPS+M algorithm: \big{(}({\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N_{0}}, (\accentset{\sim}{\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N},({\mathbf{Z}}_{\tau_{1}}^{j})_{j\leq N_{1}}\,,\dots,\,\big{)}. So when $p$ is an even number such that $p=2k$, $\mathcal{F}_{p}$ is adapted to the vector \big{(}({\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N_{0}}, $(\accentset{\sim}{\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N},\,$ \dots,\,({\mathbf{Z}}_{\tau_{k}}^{j})_{j\leq N_{k}}\big{)}. For an odd number $p=2k+1$, $\mathcal{F}_{p}$ is adapted to the vector \big{(}({\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N_{0}},(\accentset{\sim}{\mathbf{Z}}_{\tau_{0}}^{j})_{j\leq N}, \dots,\,({\mathbf{Z}}_{\tau_{k}}^{j})_{j\leq N_{k}},(\accentset{\sim}{\mathbf{Z}}_{\tau_{k}}^{j})_{j\leq N}\big{)}. For $f\in\mathscr{M}(\mathbf{E}_{\tau_{n}})$ we let $\Gamma^{\mathbf{N}}_{p,2n}(h)$ be define by

[TABLE]

and

[TABLE]

Using a telescopic argument we get

[TABLE]

with the convention for $k=0,\ \tilde{\gamma}_{-1}^{\tilde{\mathbf{N}}_{-1}}(V_{-1}Q_{0,n}(f))=\gamma_{n}(f)$.

Noticing that $\gamma_{k}^{\mathbf{N}_{k}}(1)=\tilde{\gamma}_{k-1}^{\tilde{\mathbf{N}}_{k-1}}(1)=\gamma_{k-1}^{\mathbf{N}_{k-1}}(G_{k-1})$, we can rewritte (51) as

[TABLE]

where for $k=0$, we use the convention $\gamma_{0}^{\mathbf{N}_{0}}(1)\tilde{\eta}^{\tilde{N}_{-1}}_{-1}(V_{-1}Q_{0,n}(f))=\gamma_{n}(f)$. The benefit of this decomposition is that it distinguishes the errors associated to the propagation steps and the errors associated to the selection steps. For the propagation steps, using (41) we easily get that for any $f\in\mathscr{M}(\mathbf{E}_{\tau_{k+1}})$:

[TABLE]

For the selection steps, as the resampling schemes are the same ones as for the IPS algorithm, we still have for any $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$:

[TABLE]

Thus, each selection step and propagation step is conditionally unbiased. Note that $\gamma_{k}^{\mathbf{N}_{k}}(1)$ is $\mathcal{F}_{2k-1}$-measurable and $\gamma_{k}^{\tilde{\mathbf{N}}_{k}}(1)$ is $\mathcal{F}_{2k}$-measurable, so, when the samples are generated with the IPS+M algorithm, $(\Gamma^{\mathbf{N}}_{p,2n}(h))_{p\leq 2n}$ is a $\mathcal{F}_{p}$-martingale. Therefore, $\hat{p}_{h}$ stays unbiased with the IPS+M method, because

[TABLE]

5.3.2 Almost sure convergence

Thanks to this martingale decomposition, we can use the same arguments as in the proof in the chapter 7 in [8]. The hypotheses of Theorems 7.4.2 and 7.4.3 [8, page 239 and 241] are satisfied with the IPS+M method too, which yields the following theorem:

Theorem 4**.**

For any $h\in\mathscr{M}(\mathbf{E}_{\tau_{n}})$, $\hat{p}_{h}$ converges almost surely to $p_{h}$, and, for any $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$, $\eta_{k}^{N_{k}}(f)$ converges to $\eta_{k}(f)$ almost surely, $\gamma_{k}^{N_{k}}(f)$ converges to $\gamma_{k}(f)$ almost surely.

5.3.3 A Central Limit Theorem

Theorem 5**.**

If the potential functions verify (G) and the samples are generated with the IPS+M algorithm ,then we have the following convergence in distribution:

[TABLE]

where

[TABLE]

Proof.

This proof is very similar to what is done in Chapter 9 of [8]. In order to prove that $\hat{p}_{h}$ satisfies a CLT, we begin by proving that the errors associated to the selection and propagation steps are normally distributed using Lindeberg’s theorem.

For a sequence of function $(f_{k})_{k\leq 2n}$ such that $f_{2k}$ and $f_{2k+1}$ are in $\mathscr{M}(\mathbf{E}_{\tau_{k}})$, we define the sum of errors until the $p^{th}$ selection and propagation by:

[TABLE]

For $j\in\{1,\dots N\}$ we let

[TABLE]

For $k\geq 0$, $j\in\{1,\dots,N_{k}\}$ and $i\in\{0,\dots,N^{j}\}$, we consider that the indices $j_{i}$ are ordered in such way that $j_{0}>N$ and $j_{i}<N$ when $i>0$. With such indexing $\forall s\in\{1,\dots,N\}$, $\exists j\in\{1,\dots,N_{k}\}$ and $i\in\{1,\dots,N^{j}\}$ such that $s=j_{i}$, and for such $s$ we let

[TABLE]

For $j\in\{1,\dots,N_{0}\}$, let

[TABLE]

Thus,

[TABLE]

Noting $\mathcal{P}_{k}^{N}$ a filtration adapted to the $k$ first trajectories generated in the IPS+M algorithm. Note that we have that $\mathbb{E}\left[U_{k}^{N}(f)|\mathcal{P}_{k-1}^{N}\right]=0$, and $\mathbb{E}\left[U_{k}^{N}(f)^{2}|\mathcal{P}_{k-1}^{N}\right]<\infty$, and $|U_{k}^{N}(f)|<\frac{2}{\sqrt{N}}\underset{k\leq n,\mathbf{Z}_{\tau_{k}}\in\mathbf{E}_{\tau_{k}}}{\sup}\{|f_{2k}(\mathbf{Z}_{\tau_{k}})|\wedge|f_{2k+1}(\mathbf{Z}_{\tau_{k}})|\}$, so the Lindeberg condition is clearly satisfied. Then, we have that

[TABLE]

As $\eta_{k}^{N_{k}}$ and $\tilde{\eta}_{k}^{N}$ converge almost surely to $\eta_{k}$ and $\tilde{\eta}_{k}$, $\langle\sqrt{N}M^{\mathbf{N}}_{p,2n}(f)\rangle_{n}$ converge in probability to

[TABLE]

By application of the Lindeberg’s theorem for triangular array (see for instance Theorem 4 on page 543 in [18]), we get that $\sqrt{N}M^{\mathbf{N}}_{p,2n}(f)$ converges in law to a centered Gaussian of variance $\sigma_{p}^{2}(f)$. As a corollary, if for $p\neq 2k$ we take $f_{p}=0$ and for $p=2k$ $f_{2k}=Q_{k,n}(f_{h})$, we get that

[TABLE]

and if for $p\neq 2k-1$ we take $f_{p}=0$ and for $p=2k-1$ $f_{2k-1}=V_{k-1}Q_{k,n}(f_{h})$, we get that

[TABLE]

Following from Theorem 4, \gamma_{k}^{N_{k}}(1)^{2}\ and $\tilde{\gamma}_{k}^{N}(1)^{2}$ converges almost surely to $\gamma_{k}(1)^{2}$ and $\tilde{\gamma}_{k}(1)^{2}$ , by an application of Slutsky’s Lemma, we get that $\sqrt{N}\Gamma^{\mathbf{N}}_{2N,2n}(f_{h})$ converges in law to a centered Gaussian with variance

[TABLE]

∎

5.3.4 Variance reduction

Theorem 6**.**

The variance of the original IPS can be decomposed as follows:

[TABLE]

where v_{k}(\mathbf{Z}_{\tau_{k\text{-}1}})=V_{k\text{-}1}(\mathbf{a}_{\tau_{k}}|\mathbf{Z}_{\tau_{k\text{-}1}})\big{(}1-V_{k\text{-}1}(\mathbf{a}_{\tau_{k}}|\mathbf{Z}_{\tau_{k\text{-}1}})\big{)}. Therefore we have $\sigma_{M,G}^{2}\leq\sigma_{IPS,G}^{2}$.

Proof.

[TABLE]

Temporarily using the notation $V_{k\text{-}1}(\mathbf{a}_{\tau_{k}}|\mathbf{Z}_{\tau_{k\text{-}1}})=p_{k}$, for any $f\in\mathscr{M}(\mathbf{E}_{\tau_{k}})$, we get

[TABLE]

In particular, for $f=Q_{k,n}(f_{h})$ we get

[TABLE]

Plugging (68) into the second line of (66), yields (64). ∎

6 Efficient extension of the trajectories using the Memorization method

This section presents the memorization method that was first introduced in [15]. Remember that we considered that a trajectory $\mathbf{a}_{t}$ is preponderant whenever $p_{\mathbf{a}_{t}}=\mathbb{P}(\mathbf{Z}_{t}=\mathbf{a}_{t})>0$. Assuming we know such a preponderant trajectory $\mathbf{a}_{t}$, the memorization method allows to generate a trajectory $\mathbf{Z}_{t}$ which differs from this preponderant trajectory $\mathbf{a}_{t}$.

6.1 Advantage of Memorization over a rejection algorithm

The interest of the method, compared to a rejection algorithm, is that we generate a trajectory $\mathbf{Z}_{t}\neq\mathbf{a}_{t}$ in one shot, whereas a rejection algorithm may generate several times the preponderant trajectory $\mathbf{a}_{t}$ before generating a trajectory different from $\mathbf{a}_{t}$. This is especially interesting when the probability $p_{\mathbf{a}_{t}}=\mathbb{P}(\mathbf{Z}_{t}=\mathbf{a}_{t})$ is close to 1, as, with a rejection algorithm, the average number of tries to get a trajectory different from $\mathbf{a}_{t}$ would be $\frac{1}{1-p_{\mathbf{a}_{t}}}$ which is then very high. Therefore with a rejection algorithm much computational effort would be wasted generating $\mathbf{a}_{t}$ over and over.

6.2 Remarks on improving the IPS+M

Note that the IPS+M, greatly unbalances the weights of the propagated samples. Consequently it is useless to consider an algorithm which triggers a resampling step, when the value of on the effective sample size is below a threshold. Indeed, the weights are so unbalanced that the effective sample size would always be very low and would trigger a resampling at each step.

6.3 Principle of the memorization method

6.3.1 Work with the differentiation time

The key idea of the memorization method is to consider the stopping time $\tau$ defined such that:

[TABLE]

This time $\tau$ is the time at which the trajectory $\mathbf{Z}_{t}$ differentiates itself from $\mathbf{a}_{t}$. So, to generate $\mathbf{Z}_{t}$ knowing $\tau\leq t$ is equivalent to generate $\mathbf{Z}_{t}$ knowing it differs from $\mathbf{a}_{t}$. In order to simulate a trajectory $\mathbf{Z}_{t}$ avoiding $\mathbf{a}_{t}$, one can follow these three steps:

generate $\tau$ knowing $\tau\leq t$, and set $\mathbf{Z}_{\tau^{-}}=\mathbf{a}_{\tau^{-}}$, 2. 2.

generate $Z_{\tau}$ knowing $Z_{\tau}\neq a_{\tau}$, 3. 3.

generate the rest of the trajectory normally until $t$.

These steps are not difficult to realize, except for the first one.

6.3.2 Generate $\tau$ knowing $\tau\leq t$

To achieve this first step, the authors in [15] propose to generate $\tau$ knowing $\tau\leq t$ by using a method equivalent to the inverse transform sampling method. We present hereafter the theoretical foundation for this method. We denote by $F$ the cdf of $\tau$ knowing $\tau\leq t$:

[TABLE]

and we denote by $F^{-1}$ its generalized inverse defined by

[TABLE]

We also denote by $\tilde{F}$ the function defined by

[TABLE]

where \Theta_{v}(\mathbf{a}_{v})=\big{(}(a_{s_{k}},t_{k})\big{)}_{0\leq k\leq n(\mathbf{a}_{v})}. Note that $\tilde{F}$ is discontinuous in each jump times $s_{k}$ where $K_{z_{k}^{-}}(z_{k})\neq 1$, so the inverse of $\tilde{F}$ is not necessarily defined everywhere on $[p_{\mathbf{a}_{t}},1]$. For this reason we consider $\tilde{F}^{-1}$, the generalized inverse of $\tilde{F}$ defined by

[TABLE]

$\tilde{F}^{-1}$ extends the inverse of $\tilde{F}$ constantly where it is not defined, this extension being done from the left so that $\tilde{F}^{-1}$ is right continuous.

The inverse transform sampling method consists in generating $U\sim Unif(0,1)$ and taking $F^{-1}(U)$ as a realization of $\tau\,|\,\tau\leq t$ which is a truncated random variable. The simulation of such random variables is also presented in [10]. Note that the expression of the cdf $F$ can be related to $\tilde{F}$ , indeed we have:

[TABLE]

Consequently we have that

[TABLE]

Also, as $U$ has uniform distribution on $[0,1]$, \tilde{U}=1-U\big{(}1-p_{\mathbf{a}_{t}}) is a uniform on $[p_{\mathbf{a}_{t}},1]$. Therefore, sampling with the inverse transform method is equivalent to simulating $\tilde{U}\sim Unif(p_{\mathbf{a}_{t}},1)$ and taking $\tilde{F}^{-1}(\tilde{U})$ as a realization of $\tau\,|\,\tau\leq t$.

Assuming we first generate the trajectory $\mathbf{a}_{t}$ and generate $\tilde{U}$ according to a uniform distribution on $(p_{\mathbf{a}_{t}},1)$, we now show how to evaluate $\tilde{F}^{-1}(\tilde{U})$. We consider that during the generation of $\mathbf{a}_{t}$, we computed and memorized $\mathbb{P}(\mathbf{Z}_{s_{k}^{-}}=\mathbf{a}_{s_{k}^{-}})$ and $\mathbb{P}(\mathbf{Z}_{s_{k}}=\mathbf{a}_{s_{k}})$ for each jump in the trajectory, and also for $\mathbb{P}(\mathbf{Z}_{t}=\mathbf{a}_{t})$. Then we distinguish two cases: either there exists $k\leq n(\mathbf{a}_{t})$ such that $\mathbb{P}(\mathbf{Z}_{s_{k}^{-}}=\mathbf{a}_{s_{k}^{-}})\geq\tilde{U}>\mathbb{P}(\mathbf{Z}_{s_{k}}=\mathbf{a}_{s_{k}})<$, either there exists $k\leq n(\mathbf{a}_{t})$ such that $\mathbb{P}(\mathbf{Z}_{s_{k}}=\mathbf{a}_{s_{k}})\geq\tilde{U}>\mathbb{P}(\mathbf{Z}_{s_{k}^{1}}=\mathbf{a}_{s_{k+1}^{-}})$ where we take the convention that $s_{n(\mathbf{a}_{t})+1}^{-}=t$. The first case is quite simple as by definition of $\tilde{F}^{-1}$ we get $\tilde{F}^{-1}(\tilde{U})=s_{k}$. In the second case, $\tilde{F}$ being continuous and strictly decreasing on $[s_{k},s_{k+1})$, it is inversible on this interval, and $\tilde{F}^{-1}$ corresponds to $F$’s inverse on $(\tilde{F}(s_{k+1}^{-}),\tilde{F}(s_{k})]$. So $\tilde{F}^{-1}(\tilde{U})\in[s_{k},s_{k+1})$ and $\tilde{F}(\tilde{F}^{-1}(\tilde{U}))=\tilde{U}$. Notice that

[TABLE]

So in particular, for $v=\tilde{F}^{-1}(\tilde{U})$, we have :

[TABLE]

or equivalently

[TABLE]

To determine $\tilde{F}^{-1}(\tilde{U})$ we look for the value $s$ such that the integral $\int_{0}^{s}\lambda_{a_{s_{k}}}(u)du$ is equal to $\log\left(\dfrac{\tilde{F}(s_{k})}{\tilde{U}}\right)$ by dichotomy, then we set $\tilde{F}^{-1}(\tilde{U})=s_{k}+s$.

To sum up the generation of a realization of $\tau|\tau\leq t$ we proceed as follow:

Generate $\tilde{U}\sim Unif(p_{\mathbf{a}_{t}},1)$, and set $k=0$ 2. 2.

If $\mathbb{P}(\mathbf{Z}_{s_{k}}=\mathbf{a}_{s_{k}})\geq\tilde{U}>\mathbb{P}(\mathbf{Z}_{s_{k+1}^{-}}=\mathbf{a}_{s_{k+1}^{-}})$, then we find $s\in[0,s_{k+1}-s_{k})$ such that

[TABLE]

and we set $\tau=s_{k}+s$. 3. 3.

If $\mathbb{P}(\mathbf{Z}_{s_{k+1}^{-}}=\mathbf{a}_{s_{k+1}^{-}})\geq\tilde{U}>\mathbb{P}(\mathbf{Z}_{s_{k+1}}=\mathbf{a}_{s_{k+1}})$, then $\tau=s_{k+1}$ 4. 4.

If the condition above were not satisfied, set $k=k+1$, if $k\leq n(\mathbf{a}_{s})$ repeat the steps 2 to 4

6.4 The memorization in the IPS+M algorithm

In the IPS+M algorithm, assuming we apply the memorization method on an interval $(\tau_{k\,\text{-}1},\tau_{k}]$ and for the $i^{th}$ cluster, we apply it knowing that $\mathbf{Z}_{\tau_{k\,\text{-}1}}=\mathbf{a}^{i}_{k\,\text{-}1}$. Note that this trajectory $\mathbf{a}^{i}_{k\,\text{-}1}$ is not necessarily preponderant but when we extend it into $\mathbf{a}^{i}_{k}$, the piece of trajectory $\mathbf{a}^{i}_{(\tau_{k\,\text{-}1},\tau_{k}]})$ is preponderant because :

[TABLE]

So, in the IPS we try to generate trajectories of the cluster that verify $\mathbf{Z}_{\tau_{k\,\text{-}1}}=\mathbf{a}^{i}_{k\,\text{-}1}$, but avoid the piece of trajectory $\mathbf{a}^{i}_{(\tau_{k\,\text{-}1},\tau_{k}]}$.

7 Numerical illustrations

In order to confirm our results numerically, we have applied the IPS method and the IPS+M method to two two-components system.

7.1 The heated-room system

The first system is a room heated by two heaters in passive redundancy. Heaters are programmed to maintain the temperature of the room above negative values, turning on when the temperature drops below some positive threshold and turning off when the temperature crosses a high threshold. The second heater can activate only when the first one is failed. The system fails when the temperature falls below zero.

$X_{t}$ represents the temperature of the room at time $t$. $M_{t}$ represents the status of the heaters at time $t$. Heaters can be on, off, or out-of-order, so $\mathbb{M}=\{ON,OFF,F\}^{2}$. The state of the system is $Z_{t}=(X_{t},M_{t})$.

The differential equation rule the temperature can be derived from the physics. $x_{e}$ is the exterior temperature. $\beta_{1}$ is the rate of the heat transition with the exterior. $\beta_{2}$ is the heating power of each heater. The differential equation giving the evolution of the the temperature of the room has the following form:

[TABLE]

The heaters are programmed to maintain the temperature within an interval$(x_{min},x_{max})$ where $x_{e}<0<x_{min}$. We consider that the two heaters are in passive redundancy in the sense that: when $X\leq x_{min}$ the second heater activates only if the first one is failed. When a repair of a heater occurs, if $X\leq x_{min}$ and the other heater is failed, then the heater status is set to $ON$, else the heater status is set to $OFF$. To handle the programming of the heaters, we set

$\Omega_{m}=(-\infty,x_{max})$

when all the heaters are failed $m=(F,F)$ or when at least one is activated, otherwise we set

$\Omega_{m}=(x_{min},x_{max})$

.

Due to the continuity of the temperature, the reference measure for the Kernel is $\forall B\in\mathscr{B}(E),$ $\nu_{(x,m)}(B)=\sum_{m^{+}\in\mathbb{M}\backslash\{m\}}\delta_{(x,m^{+})}(B)$. On the top boundary in $x_{max}$, heaters turn off with probability 1. On the bottom boundary in $x_{min}$, when a heater is supposed to turn on, there is a probability $\gamma=0.01$ that the heater will fail on demand. So, for instance, if z^{-}=\big{(}x_{min},(OFF,OFF)\big{)}, we have K_{z^{-}}\big{(}x_{min},(ON,OFF)\big{)}=1-\gamma, and K_{z^{-}}\big{(}x_{min},(F,ON)\big{)}=\gamma(1-\gamma), and K_{z^{-}}\big{(}x_{min},(F,F)\big{)}=\gamma^{2}.

Let $j$ be a transition from $m$ to $m^{+}$. For the spontaneous jumps that happen outside boundaries, if the transition $j$ corresponds to the failure of a heater, then: $\lambda^{j}(x,m)=0.0021+0.00015\times x$ and, if the transition corresponds to a repair, then $\lambda^{j}(x,m)=0.2\quad\mbox{when}\ M^{j}=F$. Here the system failure occurs when the temperature of the room falls below zero, so $D=\{(x,m)\in E,x<0\}$. A possible trajectory of the state of this system is plotted in figure 6. The probability of failure $p$ was estimated to $2.71\times 10^{-5}$ thanks to a massive Monte-Carlo of $10^{7}$ simulations.

7.2 Results of the simulations for the heated-room system

The results of the simulation study for the heated-room system are displayed in table 1. Here we have used the potential functions proposed in section 4. The value of $\alpha$ was set to 1.1. We have tried different values of $\alpha$ between $0.5$ and $1.5$ with a step $0.1$. The value of $\alpha=1.1$ was chosen among our trials as the one yielding the best variance reduction for the IPS method. For the IPS, IPS+M, and MC methods the variances are estimated empirically: we run the methods 100 times and we take the empirical variances of the 100 estimates. The results highlight that the IPS method is ill-suited to PDMPs, as it yields a higher variance than the MC method. Conversely, our IPS+M method performs well and has overcome the issue of the PDMP. Indeed, in the case $n=10$, it reduces the variance by a factor $2.7$ compared to the MC method, and by a factor $10$ compared to the IPS method.

The IPS+M is about 4 times slower than the IPS method, so, in terms of computational cost, the method is only $2.5$ more efficient than the IPS method on this test case. For a run of $N=10^{5}$ the IPS+M is about $2.7$ time slower than Monte-carlo method. So in terms of computational cost the IPS+M is slightly more efficient than the Monte-Carlo Method.

7.3 Remark on the SMC with Memorization

We have seen that it is possible to improve the IPS method to make it similar to the SMC method. We may, therefore, think that the IPS+M algorithm could be improved by adding adaptive optional re-sampling steps in order to get a SMC+M algorithm. In practice, however, it is not beneficial to add these adaptive optional re-sampling steps. Indeed we have noticed that, as we greatly modify the propagation process, the weights are greatly imbalanced and the effective-sample-size ends up being extremely small, which would trigger the re-sampling each time. Therefore adding adaptive optional re-sampling to the IPS+M has no effect, and in practice the IPS+M methods and the SMC+M methods are the same.

7.4 A dam system

The second system models a dam subjected to an incoming water flow. The physical variable of interest is the water level in the dam denoted by $X_{t}$. The failure of the system occurs when the water level exceeds a security threshold $x_{lim}=10$ before time $t_{f}=50$. The initial level is set to $X_{0}=0$. The water flow is characterized by the input debit $Q=10$. The dam has two evacuation valves with output debit $Q$. Each valve can be either open, close or stuck closed. So $\mathbb{M}=\{Open,Closed,Stuckclosed\}^{2}$. The valves are programmed in passive redundancy, so if the valves are in functioning order there is always one valve open and one valve closed. Though, the valve can get stuck closed and this happens at random times with exponential distribution with intensity $\lambda=0.001$. The valves are repaired with a repair rate $\mu=0.1$. When both valves are stuck closed the reservoir of the dam starts filling up according to the equation $\frac{dX_{t}}{dt}=Q/S$, where $S=10$ is the surface of the reservoir.

7.5 Results of the simulations for the dam system

The results of the simulation study for the dam system are displayed in table 2. Here we have used the potential functions:

[TABLE]

The value of $\alpha_{1}$ was set to $-0.9$ and the value of $\alpha_{2}$ was set to $-1$ (these are a priori guesses, we have not tried to use any optimization). For the IPS and the IPS+M methods the variances are estimated empirically: we run the methods 50 times and we take the empirical variances of the 50 estimates. The results are presented in table 2.

The results highlight that the IPS method is again ill-suited to PDMPs, as it yields a variance 30 times larger than the MC method. Our IPS+M method performs better than the IPS method as the variance is reduced by a factor 7. Yet on this example the IPS+M method has not overcome the issue of the PDMP, as its variance is 3.4 times larger than the variance of the Monte-carlo estimator. In terms of computational cost, on this example the IPS+M method was 3.6 times slower than the IPS, and 11.8 times slower then than the Monte-Carlo method. So the efficiency of the IPS+M is about 40 lower than the Monte-Carlo method. Clearly, the implementation of the IPS+M method requires a careful choice of the form of the potential functions and of their parameters.

8 Conclusion

This paper investigates the application of the IPS method to PDMPs. As the IPS method does not perform well when it is used on a concentrated PDMP, we introduce and analyze the IPS+M method, that is a modified version of the IPS that performs better with concentrated PDMP. The IPS+M method is similar to the IPS but has different propagation steps. Its propagation steps focus on clusters of identical particles rather then on particles individually. For each cluster a memorization method is used to get an empirical approximation of the distribution of the propagated cluster, which allows to greatly improve the accuracy of the method. We have shown that the proposed algorithm yields a strongly consistent estimation, and that this estimation satisfies a TCL. We prove that the asymptotic variance of the IPS+M estimator is always smaller than the asymptotic variance of the IPS estimator. Simulations also confirm these results, showing that the IPS+M can yield a variance reduction when the IPS cannot. In terms of computational cost, our implementations of the IPS+M method give approximately the same efficiency as the Monte-Carlo method in the examples considered in this paper, where the goal is to estimate a probability of the order of $10^{-5}$ for rather simple toy models. The numerical implementations certainly deserve more careful attention. We also believe that there are ways to improve the efficiency of the IPS+M method by finding a better class of potential functions. Another interesting improvement to the IPS+M method would be to propose an estimator of the variance. We believe that it should be possible to adapt one of the estimators proposed in [17] for the IPS method in order to get an estimator of the variance for the IPS+M estimator.