Comparison of numerical methods for the assessment of production availability of a hybrid system

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Introduction
The assessment of production availability, namely the probabilistic measure of the production regularity of a system, is a major concern in reliability studies. A recent benchmark [1] has been proposed, in order to compare different numerical methods allowing for such an evaluation in case of a small but realistic industrial situation. Such a benchmark has already been studied in [2,3]. The problem is to assess the availability, for a system of gas production, to deliver the nominal production it is meant to do, as well as other quantities, such as the annual mean number of loss of nominal production. The gas production device is what is called a hybrid system, in the sense that its time evolution is governed by two different types of dynamics: a discrete dynamic, which is related to the existence of discrete events such as failures of components and a continuous dynamic, linked to the evolution of a continuous characteristic, here the liquid level in a tank.
Such hybrid systems are typical of dynamic reliability or dynamic PSA, see [4] for instance. They may be modelled with stochastic Petri nets (see [5] with references therein) or piecewise deterministic Markov processes (PDMP), see [6] or [2]. Their numerical assessment is frequently done through Monte Carlo (MC) simulations, see [7] or [8] with references therein. It may also be done through cell to cell mapping techniques (CCMT) where the evolution of the system is approximated by a Markov chain with finite state space, see [9,10] e.g. Recently, it has been shown that finite volume (FV) methods could also be used within the same context [11]. Both methods (CCMT and FV) have already proved their capacity to be competitive with MC simulations from a computing time point of view, at least in case of small systems, see [6,10,12]. Though emerging from different scientific communities, FV methods and continuous-time CCMT (as presented in [10]) actually appear as very near: indeed, both methods start from a system of integro-differential equations fulfilled by the marginal distributions of the involved process; dividing the state space into cells, such equations are then discretized both in time and in space; such discretized equations provide what is called a FV scheme in [11] and an approximating Markov chain derived from what is called CCMT in [10]. Note, however, that the integro-differential equation used in [10] is the classical Chapman-Kolmogorov equation (as in [11]) whereas we use here another system of equations.
The aim of this paper is twofold. On one hand, we present an extension of the FV method from [11] which allows for handling such problems as the benchmark. On the other hand, in order to validate the results of the FV method, we provide confidence bands for MC simulation using a regenerative property of the system.
The paper is therefore organized as follows. We first recall the problem to be studied in Section 2 and we present the numerical results provided by regenerative Monte Carlo (RMC) simulation and the FV method. A mathematical formulation of the benchmark problem and of the different probabilistic quantities to assess is then given in Section 3 in term of a PDMP. We next recall in Section 4 the mathematical background which is needed to estimate the accuracy of RMC simulations. The FV scheme, and some of its basic properties, are presented in Section 5. Note that an important mathematical work remains to be done in order to show the convergence of the scheme towards the marginal distribution of the PDMP. Finally, we draw some future research lines in Section 6.

Description of the benchmark
For the sake of completeness, we recall here the benchmark proposed in [1] and studied in [2,3]. A gas production device is composed of two parallel production units, denoted by U 1 and U 2 , which can be up or down. Maximal production rates for the two units, respectively, are f 1;max ¼ 3200 m 3 =h and f 2;max ¼ 5500 m 3 =h. The device is required to produce gas at the nominal rate f nom ¼ 7500 m 3 =h which cannot be produced by a single unit. When both units are up, the production is then dispatched between the two of them. In order to prevent production to be under the nominal level, a reservoir R is used, with maximal capacity R ¼ 220 000 m 3 : when one or two units are down, the production is completed by taking in R the required complementary production as long as R is not empty. When one unit is down, the other one produces at its maximum level. The production rate then depends on the current volume in the reservoir: as long as R is not empty, the production is nominal. When R comes to be empty with one unit down, the production is reduced to the maximum rate of the up unit. Similarly, when both units are down, the production is nominal as long as R is nonempty and is zero after. Different cases are envisioned for the re-filling of the reservoir, which are described below.
All changes of rates of production for the units and of refilling-emptying for the reservoir are assumed to be instantaneous. Failure rates for the units are constant. Recall that for a constant failure rate l, the time to failure for a unit is exponentially distributed with the following probability density function (p.d.f.): The numerical values for the failure rates differ according to the following cases: Case 1: We first consider that the reservoir instantaneously recovers its maximal capacity R as soon as both units are up, namely at the end of repair of one unit, the other one being in up state. The respective failure rates are l 1 ¼ 1=20 000 h À1 and l 2 ¼ 1=4000 h À1 .
Case 2: We keep the same assumption as in case 1 but the failure rates are multiplied by 10.
Case 3: We keep the same assumption as in case 1 but the failure rates now depend on the production level: when a unit U i produces gas at its nominal level, namely when both unit are up (and consequently R full due to the assumption), its failure rate is the same l i as in case 2. When a unit U i produces gas at its maximum level, namely when the other unit is down, the failure rates are multiplied by 10 again (denoted by l 0 i ). Case 4: We finally remove the assumption on R instantaneously full as soon as both units are up. Also, we keep the assumption of failure rates depending on the production level, with the same numerical values as in case 3: when both units are up, if R is full, both units produce gas at their nominal level with failure rates l 1 and l 2 . If R is not full, both units are required to produce at their maximum level in order to refill R with the complementary production (rate for refilling R: f 1;max þ f 2;max À f nom ). The failure rates then are l 0 1 and l 0 2 . When unit U i is up with the other unit down, its failure rate is l 0 i too. We sum up the different failure rates in Table 1, letting l 0 i ¼ l i in the first two cases so that the same model may be used in the following for the first three cases.
Repair times are log-normally distributed, with p.d.f. f ðt i ;s i Þ ðtÞ: with the same following parameters in the four cases (note the time dimension in hours) ( Table 2). The corresponding hazard rates are then given by Case The quantities of interest for this benchmark are the following: the asymptotic availability ðA 1 Þ, which corresponds to the asymptotic ratio of time spent in nominal production, the asymptotic production availability ðPA 1 Þ, which corresponds to the asymptotic ratio between the real and the nominal productions, the asymptotic annual number (or stationary frequency) of total loss of production ðf TLP Þ, the asymptotic annual number (or stationary frequency) of loss of nominal production ðf LNP Þ.

Computational results
For each of the four envisioned cases, successive entrances in the perfect working state (where both units are up and the reservoir is full) appear as renewal times so that the system appears as a regenerative one. This allows us to compute the different quantities of interest with what we call ''regenerative'' Monte Carlo (RMC) simulations, where a sample is a cycle of the regenerative associated process (see Section 4 for details). We simulate n i.i.d. cycles in each case, with n ¼ 5 Â 10 7 in cases 1 and 2, n ¼ 10 8 in case 3 and n ¼ 10 7 in case 4. This provides us with 95% confidence bands for the different numerical results and allows us to test the results of the FV method, developed in Section 5. For each of the computational results, we have expressed the difference d between the values given by the FV method with the value given by the RMC method, in the unit s, equal to half the width of the confidence band.
We first give in Table 3 the results for the asymptotic availability A 1 .
We observe that the results given by the FV method are quite accurate. In Table 4, we give the asymptotic production availability PA 1 , and the same observation holds.
We now come to the results for the asymptotic annual numbers (or stationary frequency) of total loss of production and of nominal production. In finite horizon, such quantities are not reachable with the FV scheme. However, other quantities are accessible, which coincide with those goal quantities in long-time run. This allows to assess such stationary frequencies (see Section 3 for details).
The stationary frequency of total loss of production f TLP is provided in Table 5.
The accuracy of the FV scheme is here again quite acceptable.
We finally give the stationary frequency of loss of nominal production f LNP in Table 6.   The accuracy of the FV scheme is here not so good as for the previous computational results, but remains acceptable.

A mathematical formulation of the problem
Before developing the numerical methods used to get the previous computational results, we first specify the mathematical model for the benchmark, using notations from dynamic reliability. The time-evolution of the system is described by a PDMP ðI t ; X t Þ tX0 (see [13] or [14]). The first part I t is discrete with values in E ¼ fð1; 1Þ; ð1; 0Þ; ð0; 1Þ; ð0; 0Þg, and describes the state of both units: The second part X t stands for ''environmental condition'' according to the vocabulary of dynamic reliability. We here ; R and where, for i ¼ 1; 2, symbol X i;t stands for the duration spent by unit U i in its current state (up or down) at time t. For instance, if unit U 1 is down at time t, symbol X 1;t stands for the duration elapsed since the beginning of its repair. As for symbol X 3;t , it is the volume of gas in the reservoir at time t, with 0pX 3;t pR. Between events like failure or end of repair of one component, the discrete part I t is constant whereas the continuous part X t has a deterministic evolution (namely non-random). More precisely, assume T n and T nþ1 to be successive jumps for I t with I t ¼ Z for T n ptoT nþ1 , and X T n ¼ x 0 ¼ ðx 1;0 ; x 2;0 ; x 3;0 Þ, we then have X t ¼ g ðZÞ ðx 0 ; t À T n Þ, where g ðZÞ ðx; tÞ ¼ ðg ðZÞ 1 ðx; tÞ; g ðZÞ 2 ðx; tÞ; g ðZÞ 3 ðx; tÞÞ is given, for all ðx ¼ ðx 1 ; where v ðZÞ depends on the discrete state Z, and P ½0;R ðrÞ ¼ r if r 2 ½0; R, 0 if ro0 and R if r4R (the gas volume in the reservoir R is non-negative and cannot exceed R). Note that the function g ðZÞ satisfies the necessary semi-group property (see [14]) g ðZÞ ðx; tÞ ¼ g ðZÞ ðg ðZÞ ðx; sÞ; t À sÞ, all 0osot. The transition rate from state Z to state z at time t depends on X t ¼ x and is a function a ðZ;zÞ ðxÞ (with Z; z 2 E, Zaz, x 2 V ).
At jump time for I t from I t À ¼ Z to I t ¼ z, the environmental variable X t jumps simultaneously from X t À ¼ x to X t ¼ F ðZ;zÞ ðxÞ. This is equivalent to consider that the value of X t after the jump is random with the distribution given by the Dirac mass d F ðZ;zÞ ðxÞ ðdyÞ.   Table 9 Values of a ðZ;zÞ ðxÞ in cases 1-3 The speeds v ðZÞ and both of the functions a ðZ;zÞ ðxÞ and F ðZ;zÞ ðxÞ are given in Tables 7-12, according to the four envisioned cases, where we recall that l i and l 0 i are the constant failure rates for the units and where the hazard rates ðm i ðtÞÞ i¼1;2 are defined by (1) (see Section 2.1 for the numerical values).
The process ðI t ; X t Þ tX0 clearly is a regenerative process with successive entrances in state ðð1; 1Þ; RÞ as renewal times, where state ðð1; 1Þ; RÞ represents the perfect working state (both units are up and the reservoir is full). Mean cycle length is finite, so that a stationary distribution ðpðZ; dxÞÞ Z2E exists for the process ðI t ; X t Þ tX0 and the distribution rðtÞðÁ; dxÞ ¼ ðrðtÞðZ; dxÞÞ Z2E of ðI t ; X t Þ converges towards pðÁ; dxÞ ¼ ðpðZ; dxÞÞ Z2E . Using classical ergodic theorems, we may now express the four quantities of interest for the benchmark with p (see [15] for details on regenerative technics): Asymptotic availability: The production is nominal as long as the reservoir is non-empty (X 3;s 40). The asymptotic availability then is the asymptotic ratio between the time spent in the event fX 3;s 40g on ½0; t and t: namely: where for any event A, we set Using ergodic theorem, we get the existence of the limit (almost surely) and, for almost all o 2 O: Denoting a.s. for almost surely, we derive: because R is almost surely non-empty in state ð1; 1Þ so that pðð1; 1Þ; R 2 þ Â f0gÞ ¼ 0. Indeed, as soon as both units are up, the reservoir is either instantaneously full (cases 1-3) or, if it is empty (which may happen only in case 4), its refilling instantaneously begins.
Asymptotic production availability: Let f s be the production rate at time s. When R is non-empty ðX 3;s 40Þ, the production is nominal ðf s ¼ f nom Þ. The production is f 1;max when unit U 1 is up, unit U 2 is down and R is empty (I s ¼ ð1; 0Þ and X 3;s ¼ 0), the same for f 2;max . Otherwise (U 1 and U 2 down with R empty), the produc-tion is zero. We get which gives Stationary (annual) frequency of total loss of production: There is total loss of production when the system enters the state ðð0; 0Þ; 0Þ, namely the state where both units are down and the reservoir is empty. This occurs in two situations: either both units are down and R comes to be empty (X 3;s À 40 and X 3;s ¼ 0), or R is empty with one unit up and the other one down, and the up unit comes to fail (I s À 2 fð1; 0Þ; ð0; 1Þg and I s ¼ ð0; 0Þ). We derive ð1 fI s ¼ð0;0Þ and X 3;s À 40 and X 3;s ¼0g þ 1 fX 3;s ¼0 and I s À 2fð1;0Þ;ð0;1Þg and I s ¼ð0;0Þg Þ, ð8Þ one year is 365:25 Â 24 ¼ 8766 h. The summation is here made on the almost surely finite number of s (0ospt) such that one of the different events happens at time s. Using such an expression implies to be able to evaluate the number of times X 3;s reaches its lower bound 0 on ½0; t (with I s ¼ ð0; 0Þ). There is no special problem by MC simulation. However, it is not that clear how to express directly such a number with respect of p. We then use the following trick: in infinite horizon, the evolution of the process may be considered as stationary, so that the event ''total loss of production'' happens just as often as the event ''recover for some production''. This second type of event is easier to count: it may only happen on the recovering of one unit when R is empty and both units previously down. This means that f TLP resumes to the stationary annual frequency of transitions between ðð0; 0Þ; 0Þ and fð1; 0Þ; ð0; 1Þg. Then, we only have to count jumps for the discrete part I t of the process (while R is empty), which is now easy to achieve with p. More precisely, we have: ÂrðsÞðð0; 0Þ; dx 1 Â dx 2 Â f0gÞ ds a.s.
Annual frequency of loss of nominal production: Loss of nominal production occurs when R comes to be empty, namely: 1 fX 3;s À 40 and X 3;s ¼0g .
Similarly as for the total loss of production, considering the process as stationary, the event ''R comes to be empty'' happens with the same frequency as the event ''R comes to be non-empty''. Such an event occurs on entering the state where both units are up with previously one unit down and one unit up, and R empty. Then, f LNP resumes to the annual frequency of entrances in state ð1; 1Þ from states ðð0; 1Þ; 0Þ or ðð1; 0Þ; 0Þ. We get

RMC simulation
The problem here is to provide confidence bands for asymptotic quantities by MC simulation. A first method might be to compute the requested quantities at some large time T, namely to simulate n independent samples up to time T and then compute the associated confidence bands. We have better use here a second method, based on the regenerative property of the system. Indeed, as already noted in the previous section, ðI t ; X t Þ tX0 is a regenerative process with successive entrances in state ðð1; 1Þ; RÞ as renewal times. Also, cycles are independent with finite mean length. Usual results from renewal theory then allow us to express the goal asymptotic quantities with mean values on one single cycle. We can then perform MC simulation to evaluate such quantities with i.i.d. cycles as samples. The expressions we start with for such computations are (4,6,8,10). As for the annual frequency of loss of nominal production for instance, we might have started from (11). However, such an expression implies numerical integration of m 1 ðX 1;s Þ and of m 2 ðX 2;s Þ on some part of the cycle, which slows down the computation. Similar a remark is valid for the annual frequency of total loss of production.
We first note that each goal quantity may be written as: where c is measurable, non-negative and ''additive'' in the sense that: cððI u ; X u Þ 0ouptþs Þ ¼ cððI u ; X u Þ 0oupt Þ þ cððI u ; X u Þ touptþs Þ for all s, t40.
Let N t be such that t N t ptot N t þ1 (and consequently N t $t=Eða 1 Þ for large t), we get Beside, cððI u ; X u Þ t N t oupt Þ=t ! 0 when t ! þ1 because for large t and lim t!þ1 Y N tþ1 =N tþ1 ¼ lim n!þ1 Y n =n ¼ 0 due to EðY n Þo þ 1.
We derive that where the Y k are independent and identically distributed (i.i.d.). Using again the fact that N t $t=Eða 1 Þ for large t, the law of large numbers then implies that: Introducing the following empirical means: a natural estimator for c 1 then is c n ¼Ȳ n =ā n with c n ! EðY 1 Þ=Eða 1 Þ ¼ c 1 a.s. when n ! þ1.
In order to compute a 1 À e confidence band I e for c 1 (Pðc 1 2 I e Þ ¼ 1 À e), we introduce Z k ¼ Y k À c 1 a k for kX0 with EðZ k Þ ¼ 0 and we apply the central limit theorem (CLT) to the i.i.d. ðZ n Þ n2N . (An alternative would be to apply the CLT to ða k ; Y k Þ kX1 and then use the transformation ða; yÞ7 À!y=a, see [15]). With that aim, we need an estimator for varðZ 1 Þ with We then introduce the (unbiased) empirical variances and covariances: with S ðY Þ n ! var ðY 1 Þ a.s., S ðaÞ n ! var ða 1 Þ a.s. and S ða;Y Þ n ! covar ða 1 ; Y 1 Þ a.s. when n ! þ1, and consequently, setting S ðZÞ n ¼ S ðY Þ n À 2c n S ða;Y Þ n þ c 2 n S ðaÞ n , we have S ðZÞ n ! var ðZ 1 Þ a.s. when n ! þ1. Let w be the 1 À ðe=2Þ quantile of Nð0; 1Þ, we now derive from the CLT that for large n: SubstitutingZ n byȲ n À c 1ān and using c n ¼Ȳ n =ā n , we easily get that ½c n À s n ; c n þ s n is a 1 À e confidence band for c 1 with For the numerical results given in Section 2.2, we take e ¼ 5% and w ¼ 1:96. To get the 95% confidence band (namely c n and s n for each quantity of interest), we have to compute P n p¼1 a p , P n p¼1 a 2 p (needed for each quantity) and P n p¼1 Y p , P n p¼1 Y 2 p and P n p¼1 a p Y p for each quantity. This leads to the computation of 2 þ 4 Â 3 ¼ 14 summations in the simulations to get the four confidence bands.
The number of cycles used for the computations as well as the average cycle lengths and the computing times are given in Table 13.
We can see in such a table that the difference of models between cases 2, 3 and 4 is not very influent on the mean cycle length. As expected, a cycle in case 1 is roughly 10 times longer than in case 2 or 3. This is easily explained by the fact that the units spend most of their time in up-state with failure rates which are multiplied by 10 from the case 1 to the cases 2 or 3. Also, a cycle in case 4 is made longer by the necessity to refill the reservoir before reaching the perfect working state ðð1; 1Þ; RÞ.
The computing times are rather long but might easily be shortened improving the implementation and/or using a parallel computer. Note, however, that the computation of the confidence bands makes them clearly longer due to the necessity of computing 14 summations whereas only four summations would be required otherwise.

The FV scheme
We finally come to the main part of this paper with the presentation of the FV scheme. Such schemes are not classically used in the framework of probabilistic studies, since they have mainly been developed by engineers, in order to approximate the solutions of balance equations of the type In this equation, d is the dimension of the space variable ðx 1 ; . . . ; x d Þ, (frequently, d ¼ 2 or 3), A is a given expression of some unknown functions of the space and time variables, and the functions F i are the components of the flux of the quantity A. Then, the space domain, in which the balance equation holds, is discretized into a partition M of subdomains K 2 M called control volumes. The time is discretized using a time step dt and the above continuous equation is replaced by the following discrete one In the above equation, called the FV scheme, m K is the d-dimensional measure of K, A ðKÞ n is the approximation of A in K at the time ndt, and F n ðK; LÞ is an approximation of the flux from K to any other control volume L at the time ndt. This approximation is given as a function of the unknowns in several control volumes neighbours of K and L, which generally cannot been obtained using a mathematical procedure. On the contrary, the expression of such a function is usually based on engineering considerations, which makes higher the mathematical difficulties for a convergence proof.
We consider in this section a scheme inspired by these principles, used to directly approximate the marginal distributions of the PDMP presented in Section 3 (this is an essential difference with the MC method). The main difference between the equations satisfied by these marginal distributions and physical balance laws is that the unknown of the problem is a family of probability measures instead a family of functions. An example of the mathematical study of such a scheme in a probabilistic framework is given in [11]. Indeed, the family of probability measures ðrðtÞðZ; dxÞÞ t2R þ of the process ðI t ; X t Þ t2R þ on E Â V happens to be solution of a system of integro-differential equations [17]. As already told in the Introduction, we do not use here the classical Chapman--Kolmogorov equations as in [10,11] and we prefer use other equations which are obtained looking at the different flux of probability mass on some time interval ½t; s [17]. Denoting by B V the set of all measurable subsets of V, and for all A 2 B V , by T ðZÞ where b ðZÞ ðxÞ ¼ P z2E a ðZ;zÞ ðxÞ. The first term in the right hand of (13) corresponds to the transport of the probability mass along the characteristics u7 À!g ðZÞ ðx; uÞ on ½t; s, the second term to the exit of probability mass due to jumps out of Z on ½t; s and the third term to the arrival of probability mass due to jumps towards Z on ½t; s.
To define the FV scheme for these equations, we first define a family ðM ðZÞ Þ Z2E , such that, for all Z 2 E, M ðZÞ is a finite partition of ½0; T ðZÞ 1 Â ½0; T ðZÞ 2 Â ½0; R by polyhedral subsets, for large enough times T ðZÞ 1 and T ðZÞ 2 so that most part of the mass of rðTÞðZ; dxÞ is in ½0; T ðZÞ 1 Â ½0; T ðZÞ 2 Â ½0; R. We then select a value dt40 as time step. Following the ideas given in the introduction of this section, we intend to define an approximation P M;dt where ðp ðZ;KÞ n Þ n2N;Z2E;K2M ðZÞ is a family of real values to be constructed, recursively on n. The asymptotic distribution pðZ; dxÞ is then approximated by the FV approximation of rðTÞðZ; dxÞ for T large enough so that the distribution is nearly stabilized.
In order to construct the FV scheme (we consider in this paper an implicit version of the FV scheme given in [11]), we start from (13) divided by dt, in which symbol s is replaced by t þ dt, symbol A by K 2 M ðZÞ . We also use the relation R . We thus get  where the vector field v ðZÞ is regular enough (this is assumed in [11], but is not true on the whole domain V in the case handled in this paper), and assuming that LaK is such that K \ La; is regular enough, then v ðZÞ KL is accurately approximated bȳ where n KL ðxÞ is the unit vector, normal to K \ L and oriented from K to L, and dsðxÞ is the 2-dimensional Lebesgue measure on K \ L. Such an expression may then be used to compute v ðZÞ KL . We now turn to the approximation of term T 3 , which could simply be done by e We finally turn to the approximation e T 4 of the term T 4 , which we chose such that there is conservation of the probability mass. Using the previous approximation q ðZ;KÞ nþ1 of the probability mass escaping from fZg Â K to E Â V (see term T 3 ), we take: The scheme is now obtained by e T 1 þ e T 2 ¼ À e T 3 þ e T 4 , which gives, using the expressions given by (16) Note that it is easy to solve (24) for the particular data handled in this paper, by numbering the control volumes K 2 M ðZÞ in each state Z 2 E in such a way that v ðZÞ KL 40 implies that the number of L is greater than that of K. Therefore (24) is a simple triangular linear system with unknowns p ðZ;KÞ ðlþ1Þ . It is then possible to prove that the above algorithm converges, since it leads to a contraction in some discrete space. However, in the case considered in this paper, some faster methods (such as the BiCGSTAB algorithm) have been used.
The numerical results given in Section 2.2 have been computed using expressions (5,7,9,12) of the goal quantities substituting pðZ; dxÞ by P M;dt T ðÁ; xÞ dx for some large T. The value of T has been chosen equal to 100 000 h in cases 1 and 2, and to 50 000 h in cases 3 and 4. In each case, a grid 110 Â 110 Â 100 has been used for partition of ½0; T ðZÞ 1 Â ½0; T ðZÞ 2 Â ½0; R. We have set, in the four cases, T ð1;1Þ ¼ 10 5 h. The computing times vary from 9 mn in case 1 to less than one hour in case 4. Note that, just as for RMC simulations, a parallel computer might be used for the FV method, but it would require much more programming work. Therefore the comparison of both methods cannot be entirely based on the computing time criterion. However, it seems clear that the FV method may give similar results as MC simulation with much quicker computations. Note that for both methods, the precision may be improved either by increasing the number of cycles for MC or by adjusting the different parameters for the FV scheme.

Conclusion
This paper shows that the FV method can be successfully used to approximate the marginal distributions of stochastic processes involved in the model of hybrid systems, as an alternative to MC simulations. It leads to accurate results for an admissible computing time. Nevertheless, some mathematical work remains to be performed as the proof of the convergence of the numerical scheme. Also, some computational work must be done in order to make possible the handling of high-dimensional problems.